/*

FitModeL :  a program that Fits substitution models to the data
using Maximum Likelihood.  

Copyright (C) Stephane Guindon. jul. 2004 onward

All parts of  the source except where indicated  are distributed under
the GNU public licence.  See http://www.opensource.org for details.

*/

#include "config.h"
#include "utilities.h"
#include "models.h"
#include "eigen.h"
#include "free.h"
#include "eigenmb.h"
#include "options.h"

/*********************************************************/

int PMat_TN93(fit_double l, model *mod, fit_double **Pij)
{
  /* Analytical formula for the change probabilities
     under Tamura and Nei (1993) model.
  */

  int i,j;
  fit_double e1,e2,e3;
  fit_double a1t,a2t,bt;
  fit_double A,C,G,T,R,Y;
  fit_double kappa1,kappa2;
  
  A = mod->pi[0]; C = mod->pi[1]; G = mod->pi[2]; T = mod->pi[3];
  R = A+G;  Y = T+C;
    
  if(mod->kappa < .0) mod->kappa = 1.0e-5;
  
  if(mod->model_number < 5) { mod->lambda = 1.; }
  else if(mod->model_number == 5)
    {
      do
	{
	  mod->lambda = (Y+(R-Y)/(2.*mod->kappa))/(R-(R-Y)/(2.*mod->kappa));
	  
	  if(mod->lambda < .0)
	    {
	      mod->kappa += mod->kappa/10.;
	    }
	}while(mod->lambda < .0);
    }
  
  kappa1 = mod->kappa*2.*mod->lambda/(1.+mod->lambda);
  kappa2 = mod->kappa*2./(1.+mod->lambda);
  
  
  bt = l/(2.*(A*G*kappa1+C*T*kappa2+R*Y));
  
  a1t = kappa1;
  a2t = kappa2;
  a1t*=bt; a2t*=bt;
  
  e1 = (fit_double)exp(-a1t*R-bt*Y);
  e2 = (fit_double)exp(-a2t*Y-bt*R);
  e3 = (fit_double)exp(-bt);
  
  
  /*A->A*/Pij[0][0] = A+Y*A/R*e3+G/R*e1;
  /*A->C*/Pij[0][1] = C*(1-e3);
  /*A->G*/Pij[0][2] = G+Y*G/R*e3-G/R*e1;
  /*A->T*/Pij[0][3] = T*(1-e3);
  
  /*C->A*/Pij[1][0] = A*(1-e3);
  /*C->C*/Pij[1][1] = C+R*C/Y*e3+T/Y*e2;
  /*C->G*/Pij[1][2] = G*(1-e3);
  /*C->T*/Pij[1][3] = T+R*T/Y*e3-T/Y*e2;
  
  /*G->A*/Pij[2][0] = A+Y*A/R*e3-A/R*e1;
  /*G->C*/Pij[2][1] = C*(1-e3);
  /*G->G*/Pij[2][2] = G+Y*G/R*e3+A/R*e1;
  /*G->T*/Pij[2][3] = T*(1-e3);
  
  /*T->A*/Pij[3][0] = A*(1-e3);
  /*T->C*/Pij[3][1] = C+R*C/Y*e3-C/Y*e2;
  /*T->G*/Pij[3][2] = G*(1-e3);
  /*T->T*/Pij[3][3] = T+R*T/Y*e3+C/Y*e2;
  
  
  For(i,4) For(j,4)
    if(Pij[i][j] < 0.0) Pij[i][j] = 0.0;
  
  return 1;
}

/*********************************************************/

void dPMat_TN93(fit_double l, fit_double **dPij, model *mod, fit_double rr)
{

  /* Analytical formula for the first derivatives of the change
     probabilities under Tamura and Nei (1993) model.
     Deprecated.
  */

  fit_double kappa1,kappa2;
  fit_double de1dl,de2dl,de3dl;
  fit_double A,C,G,T,R,Y;
  fit_double Z;


  A = mod->pi[0]; C = mod->pi[1]; G = mod->pi[2]; T = mod->pi[3];
  R = A+G; Y = C+T;

  if(mod->model_number < 5) { mod->lambda = 1.; }
  else if(mod->model_number == 5)
    {
      mod->lambda = (Y+(R-Y)/(2.*mod->kappa))/(R-(R-Y)/2.*mod->kappa);
    }

  kappa1 = mod->kappa*2.*mod->lambda/(1.+mod->lambda);
  kappa2 = mod->kappa*2./(1.+mod->lambda);


  Z = 2.*A*G*kappa1+2.*C*T*kappa2+2.*R*Y;
  
  de1dl = -rr*(kappa1*R/Z + Y/Z)*(fit_double)exp(-kappa1*R/Z*l*rr-Y/Z*l*rr);
  de2dl = -rr*(kappa2*Y/Z + R/Z)*(fit_double)exp(-kappa2*Y/Z*l*rr-R/Z*l*rr);
  de3dl = -rr/Z*(fit_double)exp(-l*rr/Z);

  /*A->A*/dPij[0][0] = Y*A/R*de3dl+G/R*de1dl;
  /*A->C*/dPij[0][1] = -C*de3dl;
  /*A->G*/dPij[0][2] = Y*G/R*de3dl-G/R*de1dl;
  /*A->T*/dPij[0][3] = -T*de3dl;

  /*C->A*/dPij[1][0] = -A*de3dl;
  /*C->C*/dPij[1][1] = R*C/Y*de3dl+T/Y*de2dl;
  /*C->G*/dPij[1][2] = -G*de3dl;
  /*C->T*/dPij[1][3] = R*T/Y*de3dl-T/Y*de2dl;
 
  /*G->A*/dPij[2][0] = Y*A/R*de3dl-A/R*de1dl;
  /*G->C*/dPij[2][1] = -C*de3dl;
  /*G->G*/dPij[2][2] = Y*G/R*de3dl+A/R*de1dl;
  /*G->T*/dPij[2][3] = -T*de3dl;

  /*T->A*/dPij[3][0] = -A*de3dl;
  /*T->C*/dPij[3][1] = R*C/Y*de3dl-C/Y*de2dl;
  /*T->G*/dPij[3][2] = -G*de3dl;
  /*T->T*/dPij[3][3] = R*T/Y*de3dl+C/Y*de2dl;
}

/*********************************************************/

void d2PMat_TN93(fit_double l, fit_double **d2Pij, model *mod, fit_double rr)
{
  /* Analytical formula for the second derivatives of the change
     probabilities under Tamura and Nei (1993) model.
     Deprecated.
  */

  fit_double kappa1,kappa2;
  fit_double d2e1dl2,d2e2dl2,d2e3dl2;
  fit_double A,C,G,T,R,Y;
  fit_double Z;


  A = mod->pi[0]; C = mod->pi[2]; G = mod->pi[3]; T = mod->pi[4];
  R = A+G; Y = C+T;


  if(mod->model_number < 5) { mod->lambda = 1.; }
  else if(mod->model_number == 5)
    {
      mod->lambda = (Y+(R-Y)/(2.*mod->kappa))/(R-(R-Y)/2.*mod->kappa);
    }

  kappa1 = mod->kappa*2.*mod->lambda/(1.+mod->lambda);
  kappa2 = mod->kappa*2./(1.+mod->lambda);


  Z = 2.*A*G*kappa1+2.*C*T*kappa2+2.*R*Y;
  
  d2e1dl2 = (-rr*(kappa1*R/Z + Y/Z))*
            (-rr*(kappa1*R/Z + Y/Z))*
             (fit_double)exp(-kappa1*R/Z*l*rr-Y/Z*l*rr);
  d2e2dl2 = (-rr*(kappa2*Y/Z + R/Z))*
            (-rr*(kappa2*Y/Z + R/Z))*
             (fit_double)exp(-kappa2*Y/Z*l*rr-R/Z*l*rr);
  d2e3dl2 = (-rr/Z)*
            (-rr/Z)*
             (fit_double)exp(-l*rr/Z);

  /*A->A*/d2Pij[0][0] = Y*A/R*d2e3dl2+G/R*d2e1dl2;
  /*A->C*/d2Pij[0][1] = -C*d2e3dl2;
  /*A->G*/d2Pij[0][2] = Y*G/R*d2e3dl2-G/R*d2e1dl2;
  /*A->T*/d2Pij[0][3] = -T*d2e3dl2;

  /*C->A*/d2Pij[1][0] = -A*d2e3dl2;
  /*C->C*/d2Pij[1][1] = R*C/Y*d2e3dl2+T/Y*d2e2dl2;
  /*C->G*/d2Pij[1][2] = -G*d2e3dl2;
  /*C->T*/d2Pij[1][3] = R*T/Y*d2e3dl2-T/Y*d2e2dl2;

  /*G->A*/d2Pij[2][0] = Y*A/R*d2e3dl2-A/R*d2e1dl2;
  /*G->C*/d2Pij[2][1] = -C*d2e3dl2;
  /*G->G*/d2Pij[2][2] = Y*G/R*d2e3dl2+A/R*d2e1dl2;
  /*G->T*/d2Pij[2][3] = -T*d2e3dl2;

  /*T->A*/d2Pij[3][0] = -A*d2e3dl2;
  /*T->C*/d2Pij[3][1] = R*C/Y*d2e3dl2-C/Y*d2e2dl2;
  /*T->G*/d2Pij[3][2] = -G*d2e3dl2;
  /*T->T*/d2Pij[3][3] = R*T/Y*d2e3dl2+C/Y*d2e2dl2;

}

/*********************************************************/

int Matinv (fit_double *x, int n, int m, fit_double *space)
{
  /* Invert matrix x.  x[n*m]  ... m>=n */

   int i,j,k;
   int *irow;
   fit_double ee, t,t1,xmax;
   fit_double det;
   
   ee = 1.0E-10;
   det = 1.0;

   irow = (int *)mCalloc(n,sizeof(int));

   For (i,n)  
     {
       xmax = 0.;
       for (j=i; j<n; j++)
         if (xmax < (fit_double)fabs(x[j*m+i])) 
	   { 
	     xmax = (fit_double)fabs(x[j*m+i]); 
	     irow[i]=j; 
	   }

      det *= xmax;
      if (xmax < ee)   
	{
            printf("\n. Cannot invert the matrix of eigen vectors.\n"); fflush(NULL);
            return(0);
	}

      if (irow[i] != i) 
	{
	  For (j,m) 
	    {
	      t = x[i*m+j];
	      x[i*m+j] = x[irow[i]*m+j];
	      x[irow[i]*m+j] = t;
	    }
	}
      t = 1./x[i*m+i];
      For (j,n) 
	{
	  if (j == i) continue;
	  t1 = t*x[j*m+i];
	  For(k,m)  x[j*m+k] -= t1*x[i*m+k];
	  x[j*m+i] = -t1;
	}
      For(j,m)   x[i*m+j] *= t;
      x[i*m+i] = t;
   }                            /* i  */
   for (i=n-1; i>=0; i--) 
     {
       if (irow[i] == i) continue;
       For(j,n)  
	 {
	   t = x[j*m+i];
	   x[j*m+i] = x[j*m + irow[i]];
	   x[j*m + irow[i]] = t;
	 }
     }

   free(irow);
   return (1);
}

/*********************************************************/

int PMat(fit_double l, model *mod, fit_double **Pij, qmat *qmat_struct, arbre *tree)
{
  if(mod->model_number < 7)
  /* Use analytical formula if the substitution 
     model is nested within TN93 
  */
    return PMat_TN93(l,mod,Pij); 
  else
  /* Use numerical method otherwise */
    return PMat_Numeric(l,mod->ns,Pij,qmat_struct,tree);
}

/*********************************************************/

void dPMat(fit_double l, fit_double rr, model *mod, fit_double **dPij)
{
  /* Deprecated */
  dPMat_TN93(l,dPij,mod,rr);
}

/*********************************************************/

void d2PMat(fit_double l, fit_double rr, model *mod, fit_double **d2Pij)
{
  /* Deprecated */
  d2PMat_TN93(l,d2Pij,mod,rr);
}

/*********************************************************/

int Init_Qmat_Dayhoff(fit_double *daa, fit_double *pi)
{
  /* 
     Adapted from Ziheng Yang's PAML 

     Dayhoff's model data
     Dayhoff, M.O., Schwartz, R.M., Orcutt, B.C. (1978)
     "A model of evolutionary change in proteins."
     Dayhoff, M.O.(ed.) Atlas of Protein Sequence Structur., Vol5, Suppl3.
     National Biomedical Research Foundation, Washington DC, pp.345-352.
  */

  int i,j,naa;

  naa = 20;

  daa[ 1*20+ 0] =   27.00; daa[ 2*20+ 0] =   98.00; daa[ 2*20+ 1] =   32.00; daa[ 3*20+ 0] =  120.00;
  daa[ 3*20+ 1] =    0.00; daa[ 3*20+ 2] =  905.00; daa[ 4*20+ 0] =   36.00; daa[ 4*20+ 1] =   23.00;
  daa[ 4*20+ 2] =    0.00; daa[ 4*20+ 3] =    0.00; daa[ 5*20+ 0] =   89.00; daa[ 5*20+ 1] =  246.00;
  daa[ 5*20+ 2] =  103.00; daa[ 5*20+ 3] =  134.00; daa[ 5*20+ 4] =    0.00; daa[ 6*20+ 0] =  198.00;
  daa[ 6*20+ 1] =    1.00; daa[ 6*20+ 2] =  148.00; daa[ 6*20+ 3] = 1153.00; daa[ 6*20+ 4] =    0.00;
  daa[ 6*20+ 5] =  716.00; daa[ 7*20+ 0] =  240.00; daa[ 7*20+ 1] =    9.00; daa[ 7*20+ 2] =  139.00;
  daa[ 7*20+ 3] =  125.00; daa[ 7*20+ 4] =   11.00; daa[ 7*20+ 5] =   28.00; daa[ 7*20+ 6] =   81.00;
  daa[ 8*20+ 0] =   23.00; daa[ 8*20+ 1] =  240.00; daa[ 8*20+ 2] =  535.00; daa[ 8*20+ 3] =   86.00;
  daa[ 8*20+ 4] =   28.00; daa[ 8*20+ 5] =  606.00; daa[ 8*20+ 6] =   43.00; daa[ 8*20+ 7] =   10.00;
  daa[ 9*20+ 0] =   65.00; daa[ 9*20+ 1] =   64.00; daa[ 9*20+ 2] =   77.00; daa[ 9*20+ 3] =   24.00;
  daa[ 9*20+ 4] =   44.00; daa[ 9*20+ 5] =   18.00; daa[ 9*20+ 6] =   61.00; daa[ 9*20+ 7] =    0.00;
  daa[ 9*20+ 8] =    7.00; daa[10*20+ 0] =   41.00; daa[10*20+ 1] =   15.00; daa[10*20+ 2] =   34.00;
  daa[10*20+ 3] =    0.00; daa[10*20+ 4] =    0.00; daa[10*20+ 5] =   73.00; daa[10*20+ 6] =   11.00;
  daa[10*20+ 7] =    7.00; daa[10*20+ 8] =   44.00; daa[10*20+ 9] =  257.00; daa[11*20+ 0] =   26.00;
  daa[11*20+ 1] =  464.00; daa[11*20+ 2] =  318.00; daa[11*20+ 3] =   71.00; daa[11*20+ 4] =    0.00;
  daa[11*20+ 5] =  153.00; daa[11*20+ 6] =   83.00; daa[11*20+ 7] =   27.00; daa[11*20+ 8] =   26.00;
  daa[11*20+ 9] =   46.00; daa[11*20+10] =   18.00; daa[12*20+ 0] =   72.00; daa[12*20+ 1] =   90.00;
  daa[12*20+ 2] =    1.00; daa[12*20+ 3] =    0.00; daa[12*20+ 4] =    0.00; daa[12*20+ 5] =  114.00;
  daa[12*20+ 6] =   30.00; daa[12*20+ 7] =   17.00; daa[12*20+ 8] =    0.00; daa[12*20+ 9] =  336.00;
  daa[12*20+10] =  527.00; daa[12*20+11] =  243.00; daa[13*20+ 0] =   18.00; daa[13*20+ 1] =   14.00;
  daa[13*20+ 2] =   14.00; daa[13*20+ 3] =    0.00; daa[13*20+ 4] =    0.00; daa[13*20+ 5] =    0.00;
  daa[13*20+ 6] =    0.00; daa[13*20+ 7] =   15.00; daa[13*20+ 8] =   48.00; daa[13*20+ 9] =  196.00;
  daa[13*20+10] =  157.00; daa[13*20+11] =    0.00; daa[13*20+12] =   92.00; daa[14*20+ 0] =  250.00;
  daa[14*20+ 1] =  103.00; daa[14*20+ 2] =   42.00; daa[14*20+ 3] =   13.00; daa[14*20+ 4] =   19.00;
  daa[14*20+ 5] =  153.00; daa[14*20+ 6] =   51.00; daa[14*20+ 7] =   34.00; daa[14*20+ 8] =   94.00;
  daa[14*20+ 9] =   12.00; daa[14*20+10] =   32.00; daa[14*20+11] =   33.00; daa[14*20+12] =   17.00;
  daa[14*20+13] =   11.00; daa[15*20+ 0] =  409.00; daa[15*20+ 1] =  154.00; daa[15*20+ 2] =  495.00;
  daa[15*20+ 3] =   95.00; daa[15*20+ 4] =  161.00; daa[15*20+ 5] =   56.00; daa[15*20+ 6] =   79.00;
  daa[15*20+ 7] =  234.00; daa[15*20+ 8] =   35.00; daa[15*20+ 9] =   24.00; daa[15*20+10] =   17.00;
  daa[15*20+11] =   96.00; daa[15*20+12] =   62.00; daa[15*20+13] =   46.00; daa[15*20+14] =  245.00;
  daa[16*20+ 0] =  371.00; daa[16*20+ 1] =   26.00; daa[16*20+ 2] =  229.00; daa[16*20+ 3] =   66.00;
  daa[16*20+ 4] =   16.00; daa[16*20+ 5] =   53.00; daa[16*20+ 6] =   34.00; daa[16*20+ 7] =   30.00;
  daa[16*20+ 8] =   22.00; daa[16*20+ 9] =  192.00; daa[16*20+10] =   33.00; daa[16*20+11] =  136.00;
  daa[16*20+12] =  104.00; daa[16*20+13] =   13.00; daa[16*20+14] =   78.00; daa[16*20+15] =  550.00;
  daa[17*20+ 0] =    0.00; daa[17*20+ 1] =  201.00; daa[17*20+ 2] =   23.00; daa[17*20+ 3] =    0.00;
  daa[17*20+ 4] =    0.00; daa[17*20+ 5] =    0.00; daa[17*20+ 6] =    0.00; daa[17*20+ 7] =    0.00;
  daa[17*20+ 8] =   27.00; daa[17*20+ 9] =    0.00; daa[17*20+10] =   46.00; daa[17*20+11] =    0.00;
  daa[17*20+12] =    0.00; daa[17*20+13] =   76.00; daa[17*20+14] =    0.00; daa[17*20+15] =   75.00;
  daa[17*20+16] =    0.00; daa[18*20+ 0] =   24.00; daa[18*20+ 1] =    8.00; daa[18*20+ 2] =   95.00;
  daa[18*20+ 3] =    0.00; daa[18*20+ 4] =   96.00; daa[18*20+ 5] =    0.00; daa[18*20+ 6] =   22.00;
  daa[18*20+ 7] =    0.00; daa[18*20+ 8] =  127.00; daa[18*20+ 9] =   37.00; daa[18*20+10] =   28.00;
  daa[18*20+11] =   13.00; daa[18*20+12] =    0.00; daa[18*20+13] =  698.00; daa[18*20+14] =    0.00;
  daa[18*20+15] =   34.00; daa[18*20+16] =   42.00; daa[18*20+17] =   61.00; daa[19*20+ 0] =  208.00;
  daa[19*20+ 1] =   24.00; daa[19*20+ 2] =   15.00; daa[19*20+ 3] =   18.00; daa[19*20+ 4] =   49.00;
  daa[19*20+ 5] =   35.00; daa[19*20+ 6] =   37.00; daa[19*20+ 7] =   54.00; daa[19*20+ 8] =   44.00;
  daa[19*20+ 9] =  889.00; daa[19*20+10] =  175.00; daa[19*20+11] =   10.00; daa[19*20+12] =  258.00;
  daa[19*20+13] =   12.00; daa[19*20+14] =   48.00; daa[19*20+15] =   30.00; daa[19*20+16] =  157.00;
  daa[19*20+17] =    0.00; daa[19*20+18] =   28.00;
  
  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
  
  pi[ 0] = 0.087127; pi[ 1] = 0.040904; pi[ 2] = 0.040432; pi[ 3] = 0.046872;
  pi[ 4] = 0.033474; pi[ 5] = 0.038255; pi[ 6] = 0.049530; pi[ 7] = 0.088612;
  pi[ 8] = 0.033618; pi[ 9] = 0.036886; pi[10] = 0.085357; pi[11] = 0.080482;
  pi[12] = 0.014753; pi[13] = 0.039772; pi[14] = 0.050680; pi[15] = 0.069577;
  pi[16] = 0.058542; pi[17] = 0.010494; pi[18] = 0.029916; pi[19] = 0.064718;
 
  
    
  return 1;
}

/*********************************************************/

int Init_Qmat_DCMut(fit_double *daa, fit_double *pi)
{
  /* 
     DCMut : new implementation based on Dayhoff et al.'s raw data and amino acid mutabilities
     C. Kosiol and N. Goldman
     http://www.ebi.ac.uk/goldman-srv/dayhoff/
  */

  int i,j,naa;

  naa = 20;

  daa[ 1*20+ 0] =   26.78280; daa[ 2*20+ 0] =   98.44740; daa[ 2*20+ 1] =   32.70590; daa[ 3*20+ 0] =  119.98050; 
  daa[ 3*20+ 1] =    0.00000; daa[ 3*20+ 2] =  893.15150; daa[ 4*20+ 0] =   36.00160; daa[ 4*20+ 1] =   23.23740; 
  daa[ 4*20+ 2] =    0.00000; daa[ 4*20+ 3] =    0.00000; daa[ 5*20+ 0] =   88.77530; daa[ 5*20+ 1] =  243.99390; 
  daa[ 5*20+ 2] =  102.85090; daa[ 5*20+ 3] =  134.85510; daa[ 5*20+ 4] =    0.00000; daa[ 6*20+ 0] =  196.11670; 
  daa[ 6*20+ 1] =    0.00000; daa[ 6*20+ 2] =  149.34090; daa[ 6*20+ 3] = 1138.86590; daa[ 6*20+ 4] =    0.00000; 
  daa[ 6*20+ 5] =  708.60220; daa[ 7*20+ 0] =  238.61110; daa[ 7*20+ 1] =    8.77910; daa[ 7*20+ 2] =  138.53520; 
  daa[ 7*20+ 3] =  124.09810; daa[ 7*20+ 4] =   10.72780; daa[ 7*20+ 5] =   28.15810; daa[ 7*20+ 6] =   81.19070; 
  daa[ 8*20+ 0] =   22.81160; daa[ 8*20+ 1] =  238.31480; daa[ 8*20+ 2] =  529.00240; daa[ 8*20+ 3] =   86.82410; 
  daa[ 8*20+ 4] =   28.27290; daa[ 8*20+ 5] =  601.16130; daa[ 8*20+ 6] =   43.94690; daa[ 8*20+ 7] =   10.68020; 
  daa[ 9*20+ 0] =   65.34160; daa[ 9*20+ 1] =   63.26290; daa[ 9*20+ 2] =   76.80240; daa[ 9*20+ 3] =   23.92480; 
  daa[ 9*20+ 4] =   43.80740; daa[ 9*20+ 5] =   18.03930; daa[ 9*20+ 6] =   60.95260; daa[ 9*20+ 7] =    0.00000; 
  daa[ 9*20+ 8] =    7.69810; daa[10*20+ 0] =   40.64310; daa[10*20+ 1] =   15.49240; daa[10*20+ 2] =   34.11130; 
  daa[10*20+ 3] =    0.00000; daa[10*20+ 4] =    0.00000; daa[10*20+ 5] =   73.07720; daa[10*20+ 6] =   11.28800; 
  daa[10*20+ 7] =    7.15140; daa[10*20+ 8] =   44.35040; daa[10*20+ 9] =  255.66850; daa[11*20+ 0] =   25.86350; 
  daa[11*20+ 1] =  461.01240; daa[11*20+ 2] =  314.83710; daa[11*20+ 3] =   71.69130; daa[11*20+ 4] =    0.00000; 
  daa[11*20+ 5] =  151.90780; daa[11*20+ 6] =   83.00780; daa[11*20+ 7] =   26.76830; daa[11*20+ 8] =   27.04750; 
  daa[11*20+ 9] =   46.08570; daa[11*20+10] =   18.06290; daa[12*20+ 0] =   71.78400; daa[12*20+ 1] =   89.63210; 
  daa[12*20+ 2] =    0.00000; daa[12*20+ 3] =    0.00000; daa[12*20+ 4] =    0.00000; daa[12*20+ 5] =  112.74990; 
  daa[12*20+ 6] =   30.48030; daa[12*20+ 7] =   17.03720; daa[12*20+ 8] =    0.00000; daa[12*20+ 9] =  333.27320; 
  daa[12*20+10] =  523.01150; daa[12*20+11] =  241.17390; daa[13*20+ 0] =   18.36410; daa[13*20+ 1] =   13.69060; 
  daa[13*20+ 2] =   13.85030; daa[13*20+ 3] =    0.00000; daa[13*20+ 4] =    0.00000; daa[13*20+ 5] =    0.00000; 
  daa[13*20+ 6] =    0.00000; daa[13*20+ 7] =   15.34780; daa[13*20+ 8] =   47.59270; daa[13*20+ 9] =  195.19510; 
  daa[13*20+10] =  156.51600; daa[13*20+11] =    0.00000; daa[13*20+12] =   92.18600; daa[14*20+ 0] =  248.59200; 
  daa[14*20+ 1] =  102.83130; daa[14*20+ 2] =   41.92440; daa[14*20+ 3] =   13.39400; daa[14*20+ 4] =   18.75500; 
  daa[14*20+ 5] =  152.61880; daa[14*20+ 6] =   50.70030; daa[14*20+ 7] =   34.71530; daa[14*20+ 8] =   93.37090; 
  daa[14*20+ 9] =   11.91520; daa[14*20+10] =   31.62580; daa[14*20+11] =   33.54190; daa[14*20+12] =   17.02050; 
  daa[14*20+13] =   11.05060; daa[15*20+ 0] =  405.18700; daa[15*20+ 1] =  153.15900; daa[15*20+ 2] =  488.58920; 
  daa[15*20+ 3] =   95.60970; daa[15*20+ 4] =  159.83560; daa[15*20+ 5] =   56.18280; daa[15*20+ 6] =   79.39990; 
  daa[15*20+ 7] =  232.22430; daa[15*20+ 8] =   35.36430; daa[15*20+ 9] =   24.79550; daa[15*20+10] =   17.14320; 
  daa[15*20+11] =   95.45570; daa[15*20+12] =   61.99510; daa[15*20+13] =   45.99010; daa[15*20+14] =  242.72020; 
  daa[16*20+ 0] =  368.03650; daa[16*20+ 1] =   26.57450; daa[16*20+ 2] =  227.16970; daa[16*20+ 3] =   66.09300; 
  daa[16*20+ 4] =   16.23660; daa[16*20+ 5] =   52.56510; daa[16*20+ 6] =   34.01560; daa[16*20+ 7] =   30.66620; 
  daa[16*20+ 8] =   22.63330; daa[16*20+ 9] =  190.07390; daa[16*20+10] =   33.10900; daa[16*20+11] =  135.05990; 
  daa[16*20+12] =  103.15340; daa[16*20+13] =   13.66550; daa[16*20+14] =   78.28570; daa[16*20+15] =  543.66740; 
  daa[17*20+ 0] =    0.00000; daa[17*20+ 1] =  200.13750; daa[17*20+ 2] =   22.49680; daa[17*20+ 3] =    0.00000; 
  daa[17*20+ 4] =    0.00000; daa[17*20+ 5] =    0.00000; daa[17*20+ 6] =    0.00000; daa[17*20+ 7] =    0.00000; 
  daa[17*20+ 8] =   27.05640; daa[17*20+ 9] =    0.00000; daa[17*20+10] =   46.17760; daa[17*20+11] =    0.00000; 
  daa[17*20+12] =    0.00000; daa[17*20+13] =   76.23540; daa[17*20+14] =    0.00000; daa[17*20+15] =   74.08190; 
  daa[17*20+16] =    0.00000; daa[18*20+ 0] =   24.41390; daa[18*20+ 1] =    7.80120; daa[18*20+ 2] =   94.69400; 
  daa[18*20+ 3] =    0.00000; daa[18*20+ 4] =   95.31640; daa[18*20+ 5] =    0.00000; daa[18*20+ 6] =   21.47170; 
  daa[18*20+ 7] =    0.00000; daa[18*20+ 8] =  126.54000; daa[18*20+ 9] =   37.48340; daa[18*20+10] =   28.65720; 
  daa[18*20+11] =   13.21420; daa[18*20+12] =    0.00000; daa[18*20+13] =  695.26290; daa[18*20+14] =    0.00000; 
  daa[18*20+15] =   33.62890; daa[18*20+16] =   41.78390; daa[18*20+17] =   60.80700; daa[19*20+ 0] =  205.95640; 
  daa[19*20+ 1] =   24.03680; daa[19*20+ 2] =   15.80670; daa[19*20+ 3] =   17.83160; daa[19*20+ 4] =   48.46780; 
  daa[19*20+ 5] =   34.69830; daa[19*20+ 6] =   36.72500; daa[19*20+ 7] =   53.81650; daa[19*20+ 8] =   43.87150; 
  daa[19*20+ 9] =  881.00380; daa[19*20+10] =  174.51560; daa[19*20+11] =   10.38500; daa[19*20+12] =  256.59550; 
  daa[19*20+13] =   12.36060; daa[19*20+14] =   48.50260; daa[19*20+15] =   30.38360; daa[19*20+16] =  156.19970; 
  daa[19*20+17] =    0.00000; daa[19*20+18] =   27.93790; 
  


  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];

  pi[ 0] = 0.087127; pi[ 1] = 0.040904; pi[ 2] = 0.040432; pi[ 3] = 0.046872; 
  pi[ 4] = 0.033474; pi[ 5] = 0.038255; pi[ 6] = 0.049530; pi[ 7] = 0.088612;
  pi[ 8] = 0.033619; pi[ 9] = 0.036886; pi[10] = 0.085357; pi[11] = 0.080481;
  pi[12] = 0.014753; pi[13] = 0.039772; pi[14] = 0.050680; pi[15] = 0.069577;
  pi[16] = 0.058542; pi[17] = 0.010494; pi[18] = 0.029916; pi[19] = 0.064718; 

  return 1;


}

/*********************************************************/

int Init_Qmat_WAG(fit_double *daa, fit_double *pi)
{
  int i,j,naa;

  /* 
     WAG's model data
     Simon Whelan and Nick Goldman
     'A general empirical model of protein evolution derived from multiple
     protein families using a maximum-likelihood approach' 
     MBE (2001) 18:691-699
  */

  naa = 20;

  daa[ 1*20+ 0] =  55.15710; daa[ 2*20+ 0] =  50.98480; daa[ 2*20+ 1] =  63.53460; 
  daa[ 3*20+ 0] =  73.89980; daa[ 3*20+ 1] =  14.73040; daa[ 3*20+ 2] = 542.94200; 
  daa[ 4*20+ 0] = 102.70400; daa[ 4*20+ 1] =  52.81910; daa[ 4*20+ 2] =  26.52560; 
  daa[ 4*20+ 3] =   3.02949; daa[ 5*20+ 0] =  90.85980; daa[ 5*20+ 1] = 303.55000; 
  daa[ 5*20+ 2] = 154.36400; daa[ 5*20+ 3] =  61.67830; daa[ 5*20+ 4] =   9.88179; 
  daa[ 6*20+ 0] = 158.28500; daa[ 6*20+ 1] =  43.91570; daa[ 6*20+ 2] =  94.71980; 
  daa[ 6*20+ 3] = 617.41600; daa[ 6*20+ 4] =   2.13520; daa[ 6*20+ 5] = 546.94700; 
  daa[ 7*20+ 0] = 141.67200; daa[ 7*20+ 1] =  58.46650; daa[ 7*20+ 2] = 112.55600; 
  daa[ 7*20+ 3] =  86.55840; daa[ 7*20+ 4] =  30.66740; daa[ 7*20+ 5] =  33.00520; 
  daa[ 7*20+ 6] =  56.77170; daa[ 8*20+ 0] =  31.69540; daa[ 8*20+ 1] = 213.71500; 
  daa[ 8*20+ 2] = 395.62900; daa[ 8*20+ 3] =  93.06760; daa[ 8*20+ 4] =  24.89720; 
  daa[ 8*20+ 5] = 429.41100; daa[ 8*20+ 6] =  57.00250; daa[ 8*20+ 7] =  24.94100; 
  daa[ 9*20+ 0] =  19.33350; daa[ 9*20+ 1] =  18.69790; daa[ 9*20+ 2] =  55.42360; 
  daa[ 9*20+ 3] =   3.94370; daa[ 9*20+ 4] =  17.01350; daa[ 9*20+ 5] =  11.39170; 
  daa[ 9*20+ 6] =  12.73950; daa[ 9*20+ 7] =   3.04501; daa[ 9*20+ 8] =  13.81900; 
  daa[10*20+ 0] =  39.79150; daa[10*20+ 1] =  49.76710; daa[10*20+ 2] =  13.15280; 
  daa[10*20+ 3] =   8.48047; daa[10*20+ 4] =  38.42870; daa[10*20+ 5] =  86.94890; 
  daa[10*20+ 6] =  15.42630; daa[10*20+ 7] =   6.13037; daa[10*20+ 8] =  49.94620; 
  daa[10*20+ 9] = 317.09700; daa[11*20+ 0] =  90.62650; daa[11*20+ 1] = 535.14200; 
  daa[11*20+ 2] = 301.20100; daa[11*20+ 3] =  47.98550; daa[11*20+ 4] =   7.40339; 
  daa[11*20+ 5] = 389.49000; daa[11*20+ 6] = 258.44300; daa[11*20+ 7] =  37.35580; 
  daa[11*20+ 8] =  89.04320; daa[11*20+ 9] =  32.38320; daa[11*20+10] =  25.75550; 
  daa[12*20+ 0] =  89.34960; daa[12*20+ 1] =  68.31620; daa[12*20+ 2] =  19.82210; 
  daa[12*20+ 3] =  10.37540; daa[12*20+ 4] =  39.04820; daa[12*20+ 5] = 154.52600; 
  daa[12*20+ 6] =  31.51240; daa[12*20+ 7] =  17.41000; daa[12*20+ 8] =  40.41410; 
  daa[12*20+ 9] = 425.74600; daa[12*20+10] = 485.40200; daa[12*20+11] =  93.42760; 
  daa[13*20+ 0] =  21.04940; daa[13*20+ 1] =  10.27110; daa[13*20+ 2] =   9.61621; 
  daa[13*20+ 3] =   4.67304; daa[13*20+ 4] =  39.80200; daa[13*20+ 5] =   9.99208; 
  daa[13*20+ 6] =   8.11339; daa[13*20+ 7] =   4.99310; daa[13*20+ 8] =  67.93710; 
  daa[13*20+ 9] = 105.94700; daa[13*20+10] = 211.51700; daa[13*20+11] =   8.88360; 
  daa[13*20+12] = 119.06300; daa[14*20+ 0] = 143.85500; daa[14*20+ 1] =  67.94890; 
  daa[14*20+ 2] =  19.50810; daa[14*20+ 3] =  42.39840; daa[14*20+ 4] =  10.94040; 
  daa[14*20+ 5] =  93.33720; daa[14*20+ 6] =  68.23550; daa[14*20+ 7] =  24.35700; 
  daa[14*20+ 8] =  69.61980; daa[14*20+ 9] =   9.99288; daa[14*20+10] =  41.58440; 
  daa[14*20+11] =  55.68960; daa[14*20+12] =  17.13290; daa[14*20+13] =  16.14440; 
  daa[15*20+ 0] = 337.07900; daa[15*20+ 1] = 122.41900; daa[15*20+ 2] = 397.42300; 
  daa[15*20+ 3] = 107.17600; daa[15*20+ 4] = 140.76600; daa[15*20+ 5] = 102.88700; 
  daa[15*20+ 6] =  70.49390; daa[15*20+ 7] = 134.18200; daa[15*20+ 8] =  74.01690; 
  daa[15*20+ 9] =  31.94400; daa[15*20+10] =  34.47390; daa[15*20+11] =  96.71300; 
  daa[15*20+12] =  49.39050; daa[15*20+13] =  54.59310; daa[15*20+14] = 161.32800; 
  daa[16*20+ 0] = 212.11100; daa[16*20+ 1] =  55.44130; daa[16*20+ 2] = 203.00600; 
  daa[16*20+ 3] =  37.48660; daa[16*20+ 4] =  51.29840; daa[16*20+ 5] =  85.79280; 
  daa[16*20+ 6] =  82.27650; daa[16*20+ 7] =  22.58330; daa[16*20+ 8] =  47.33070; 
  daa[16*20+ 9] = 145.81600; daa[16*20+10] =  32.66220; daa[16*20+11] = 138.69800; 
  daa[16*20+12] = 151.61200; daa[16*20+13] =  17.19030; daa[16*20+14] =  79.53840; 
  daa[16*20+15] = 437.80200; daa[17*20+ 0] =  11.31330; daa[17*20+ 1] = 116.39200; 
  daa[17*20+ 2] =   7.19167; daa[17*20+ 3] =  12.97670; daa[17*20+ 4] =  71.70700; 
  daa[17*20+ 5] =  21.57370; daa[17*20+ 6] =  15.65570; daa[17*20+ 7] =  33.69830; 
  daa[17*20+ 8] =  26.25690; daa[17*20+ 9] =  21.24830; daa[17*20+10] =  66.53090; 
  daa[17*20+11] =  13.75050; daa[17*20+12] =  51.57060; daa[17*20+13] = 152.96400; 
  daa[17*20+14] =  13.94050; daa[17*20+15] =  52.37420; daa[17*20+16] =  11.08640; 
  daa[18*20+ 0] =  24.07350; daa[18*20+ 1] =  38.15330; daa[18*20+ 2] = 108.60000; 
  daa[18*20+ 3] =  32.57110; daa[18*20+ 4] =  54.38330; daa[18*20+ 5] =  22.77100; 
  daa[18*20+ 6] =  19.63030; daa[18*20+ 7] =  10.36040; daa[18*20+ 8] = 387.34400; 
  daa[18*20+ 9] =  42.01700; daa[18*20+10] =  39.86180; daa[18*20+11] =  13.32640; 
  daa[18*20+12] =  42.84370; daa[18*20+13] = 645.42800; daa[18*20+14] =  21.60460; 
  daa[18*20+15] =  78.69930; daa[18*20+16] =  29.11480; daa[18*20+17] = 248.53900; 
  daa[19*20+ 0] = 200.60100; daa[19*20+ 1] =  25.18490; daa[19*20+ 2] =  19.62460; 
  daa[19*20+ 3] =  15.23350; daa[19*20+ 4] = 100.21400; daa[19*20+ 5] =  30.12810; 
  daa[19*20+ 6] =  58.87310; daa[19*20+ 7] =  18.72470; daa[19*20+ 8] =  11.83580; 
  daa[19*20+ 9] = 782.13000; daa[19*20+10] = 180.03400; daa[19*20+11] =  30.54340; 
  daa[19*20+12] = 205.84500; daa[19*20+13] =  64.98920; daa[19*20+14] =  31.48870; 
  daa[19*20+15] =  23.27390; daa[19*20+16] = 138.82300; daa[19*20+17] =  36.53690; 
  daa[19*20+18] =  31.47300; 

  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];

  pi[0] = 0.0866279; pi[1] =  0.043972; pi[2] =  0.0390894; pi[3] =  0.0570451;
  pi[4] =  0.0193078; pi[5] =  0.0367281; pi[6] =  0.0580589; pi[7] =  0.0832518;
  pi[8] =  0.0244313; pi[9] =  0.048466; pi[10] =  0.086209; pi[11] = 0.0620286;
  pi[12] = 0.0195027; pi[13] =  0.0384319; pi[14] =  0.0457631; pi[15] = 0.0695179;
  pi[16] =  0.0610127; pi[17] =  0.0143859; pi[18] =  0.0352742; pi[19] =  0.0708956;

  return 1;
}


/*********************************************************/

int Init_Qmat_JTT(fit_double *daa, fit_double *pi)
{
  int i,j,naa;

  /* 
     Adapted from Ziheng Yang's PAML 
	 
     JTT's model data
     D.T.Jones, W.R.Taylor and J.M.Thornton
     "The rapid generation of mutation data matrices from protein sequences"
     CABIOS  vol.8 no.3 1992 pp275-282
  */


  naa = 20;

  daa[ 1*20+ 0] =   58.00; daa[ 2*20+ 0] =   54.00; daa[ 2*20+ 1] =   45.00; daa[ 3*20+ 0] =   81.00;
  daa[ 3*20+ 1] =   16.00; daa[ 3*20+ 2] =  528.00; daa[ 4*20+ 0] =   56.00; daa[ 4*20+ 1] =  113.00;
  daa[ 4*20+ 2] =   34.00; daa[ 4*20+ 3] =   10.00; daa[ 5*20+ 0] =   57.00; daa[ 5*20+ 1] =  310.00;
  daa[ 5*20+ 2] =   86.00; daa[ 5*20+ 3] =   49.00; daa[ 5*20+ 4] =    9.00; daa[ 6*20+ 0] =  105.00;
  daa[ 6*20+ 1] =   29.00; daa[ 6*20+ 2] =   58.00; daa[ 6*20+ 3] =  767.00; daa[ 6*20+ 4] =    5.00;
  daa[ 6*20+ 5] =  323.00; daa[ 7*20+ 0] =  179.00; daa[ 7*20+ 1] =  137.00; daa[ 7*20+ 2] =   81.00;
  daa[ 7*20+ 3] =  130.00; daa[ 7*20+ 4] =   59.00; daa[ 7*20+ 5] =   26.00; daa[ 7*20+ 6] =  119.00;
  daa[ 8*20+ 0] =   27.00; daa[ 8*20+ 1] =  328.00; daa[ 8*20+ 2] =  391.00; daa[ 8*20+ 3] =  112.00;
  daa[ 8*20+ 4] =   69.00; daa[ 8*20+ 5] =  597.00; daa[ 8*20+ 6] =   26.00; daa[ 8*20+ 7] =   23.00;
  daa[ 9*20+ 0] =   36.00; daa[ 9*20+ 1] =   22.00; daa[ 9*20+ 2] =   47.00; daa[ 9*20+ 3] =   11.00;
  daa[ 9*20+ 4] =   17.00; daa[ 9*20+ 5] =    9.00; daa[ 9*20+ 6] =   12.00; daa[ 9*20+ 7] =    6.00;
  daa[ 9*20+ 8] =   16.00; daa[10*20+ 0] =   30.00; daa[10*20+ 1] =   38.00; daa[10*20+ 2] =   12.00;
  daa[10*20+ 3] =    7.00; daa[10*20+ 4] =   23.00; daa[10*20+ 5] =   72.00; daa[10*20+ 6] =    9.00;
  daa[10*20+ 7] =    6.00; daa[10*20+ 8] =   56.00; daa[10*20+ 9] =  229.00; daa[11*20+ 0] =   35.00;
  daa[11*20+ 1] =  646.00; daa[11*20+ 2] =  263.00; daa[11*20+ 3] =   26.00; daa[11*20+ 4] =    7.00;
  daa[11*20+ 5] =  292.00; daa[11*20+ 6] =  181.00; daa[11*20+ 7] =   27.00; daa[11*20+ 8] =   45.00;
  daa[11*20+ 9] =   21.00; daa[11*20+10] =   14.00; daa[12*20+ 0] =   54.00; daa[12*20+ 1] =   44.00;
  daa[12*20+ 2] =   30.00; daa[12*20+ 3] =   15.00; daa[12*20+ 4] =   31.00; daa[12*20+ 5] =   43.00;
  daa[12*20+ 6] =   18.00; daa[12*20+ 7] =   14.00; daa[12*20+ 8] =   33.00; daa[12*20+ 9] =  479.00;
  daa[12*20+10] =  388.00; daa[12*20+11] =   65.00; daa[13*20+ 0] =   15.00; daa[13*20+ 1] =    5.00;
  daa[13*20+ 2] =   10.00; daa[13*20+ 3] =    4.00; daa[13*20+ 4] =   78.00; daa[13*20+ 5] =    4.00;
  daa[13*20+ 6] =    5.00; daa[13*20+ 7] =    5.00; daa[13*20+ 8] =   40.00; daa[13*20+ 9] =   89.00;
  daa[13*20+10] =  248.00; daa[13*20+11] =    4.00; daa[13*20+12] =   43.00; daa[14*20+ 0] =  194.00;
  daa[14*20+ 1] =   74.00; daa[14*20+ 2] =   15.00; daa[14*20+ 3] =   15.00; daa[14*20+ 4] =   14.00;
  daa[14*20+ 5] =  164.00; daa[14*20+ 6] =   18.00; daa[14*20+ 7] =   24.00; daa[14*20+ 8] =  115.00;
  daa[14*20+ 9] =   10.00; daa[14*20+10] =  102.00; daa[14*20+11] =   21.00; daa[14*20+12] =   16.00;
  daa[14*20+13] =   17.00; daa[15*20+ 0] =  378.00; daa[15*20+ 1] =  101.00; daa[15*20+ 2] =  503.00;
  daa[15*20+ 3] =   59.00; daa[15*20+ 4] =  223.00; daa[15*20+ 5] =   53.00; daa[15*20+ 6] =   30.00;
  daa[15*20+ 7] =  201.00; daa[15*20+ 8] =   73.00; daa[15*20+ 9] =   40.00; daa[15*20+10] =   59.00;
  daa[15*20+11] =   47.00; daa[15*20+12] =   29.00; daa[15*20+13] =   92.00; daa[15*20+14] =  285.00;
  daa[16*20+ 0] =  475.00; daa[16*20+ 1] =   64.00; daa[16*20+ 2] =  232.00; daa[16*20+ 3] =   38.00;
  daa[16*20+ 4] =   42.00; daa[16*20+ 5] =   51.00; daa[16*20+ 6] =   32.00; daa[16*20+ 7] =   33.00;
  daa[16*20+ 8] =   46.00; daa[16*20+ 9] =  245.00; daa[16*20+10] =   25.00; daa[16*20+11] =  103.00;
  daa[16*20+12] =  226.00; daa[16*20+13] =   12.00; daa[16*20+14] =  118.00; daa[16*20+15] =  477.00;
  daa[17*20+ 0] =    9.00; daa[17*20+ 1] =  126.00; daa[17*20+ 2] =    8.00; daa[17*20+ 3] =    4.00;
  daa[17*20+ 4] =  115.00; daa[17*20+ 5] =   18.00; daa[17*20+ 6] =   10.00; daa[17*20+ 7] =   55.00;
  daa[17*20+ 8] =    8.00; daa[17*20+ 9] =    9.00; daa[17*20+10] =   52.00; daa[17*20+11] =   10.00;
  daa[17*20+12] =   24.00; daa[17*20+13] =   53.00; daa[17*20+14] =    6.00; daa[17*20+15] =   35.00;
  daa[17*20+16] =   12.00; daa[18*20+ 0] =   11.00; daa[18*20+ 1] =   20.00; daa[18*20+ 2] =   70.00;
  daa[18*20+ 3] =   46.00; daa[18*20+ 4] =  209.00; daa[18*20+ 5] =   24.00; daa[18*20+ 6] =    7.00;
  daa[18*20+ 7] =    8.00; daa[18*20+ 8] =  573.00; daa[18*20+ 9] =   32.00; daa[18*20+10] =   24.00;
  daa[18*20+11] =    8.00; daa[18*20+12] =   18.00; daa[18*20+13] =  536.00; daa[18*20+14] =   10.00;
  daa[18*20+15] =   63.00; daa[18*20+16] =   21.00; daa[18*20+17] =   71.00; daa[19*20+ 0] =  298.00;
  daa[19*20+ 1] =   17.00; daa[19*20+ 2] =   16.00; daa[19*20+ 3] =   31.00; daa[19*20+ 4] =   62.00;
  daa[19*20+ 5] =   20.00; daa[19*20+ 6] =   45.00; daa[19*20+ 7] =   47.00; daa[19*20+ 8] =   11.00;
  daa[19*20+ 9] =  961.00; daa[19*20+10] =  180.00; daa[19*20+11] =   14.00; daa[19*20+12] =  323.00;
  daa[19*20+13] =   62.00; daa[19*20+14] =   23.00; daa[19*20+15] =   38.00; daa[19*20+16] =  112.00;
  daa[19*20+17] =   25.00; daa[19*20+18] =   16.00;
  
  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];

  pi[ 0] = 0.076748; pi[ 1] = 0.051691; pi[ 2] = 0.042645; pi[ 3] = 0.051544;
  pi[ 4] = 0.019803; pi[ 5] = 0.040752; pi[ 6] = 0.061830; pi[ 7] = 0.073152;
  pi[ 8] = 0.022944; pi[ 9] = 0.053761; pi[10] = 0.091904; pi[11] = 0.058676;
  pi[12] = 0.023826; pi[13] = 0.040126; pi[14] = 0.050901; pi[15] = 0.068765;
  pi[16] = 0.058565; pi[17] = 0.014261; pi[18] = 0.032102; pi[19] = 0.066005;
 
 return 1;
}

/*********************************************************/

int Init_Qmat_MtREV(fit_double *daa, fit_double *pi)
{
  /*
    Adapted from Ziheng Yang's PAML 
    
    Adachi, J. and Hasegawa, M. (1996) MOLPHY version 2.3: programs for
    molecular phylogenetics based on maximum likelihood.  Computer Science
    Monographs of Institute of Statistical Mathematics 28:1-150.
  */

  int i,j,naa;

  naa = 20;


  daa[ 1*20+ 0] =   23.18; daa[ 2*20+ 0] =   26.95; daa[ 2*20+ 1] =   13.24; daa[ 3*20+ 0] =   17.67;
  daa[ 3*20+ 1] =    1.90; daa[ 3*20+ 2] =  794.38; daa[ 4*20+ 0] =   59.93; daa[ 4*20+ 1] =  103.33;
  daa[ 4*20+ 2] =   58.94; daa[ 4*20+ 3] =    1.90; daa[ 5*20+ 0] =    1.90; daa[ 5*20+ 1] =  220.99;
  daa[ 5*20+ 2] =  173.56; daa[ 5*20+ 3] =   55.28; daa[ 5*20+ 4] =   75.24; daa[ 6*20+ 0] =    9.77;
  daa[ 6*20+ 1] =    1.90; daa[ 6*20+ 2] =   63.05; daa[ 6*20+ 3] =  583.55; daa[ 6*20+ 4] =    1.90;
  daa[ 6*20+ 5] =  313.56; daa[ 7*20+ 0] =  120.71; daa[ 7*20+ 1] =   23.03; daa[ 7*20+ 2] =   53.30;
  daa[ 7*20+ 3] =   56.77; daa[ 7*20+ 4] =   30.71; daa[ 7*20+ 5] =    6.75; daa[ 7*20+ 6] =   28.28;
  daa[ 8*20+ 0] =   13.90; daa[ 8*20+ 1] =  165.23; daa[ 8*20+ 2] =  496.13; daa[ 8*20+ 3] =  113.99;
  daa[ 8*20+ 4] =  141.49; daa[ 8*20+ 5] =  582.40; daa[ 8*20+ 6] =   49.12; daa[ 8*20+ 7] =    1.90;
  daa[ 9*20+ 0] =   96.49; daa[ 9*20+ 1] =    1.90; daa[ 9*20+ 2] =   27.10; daa[ 9*20+ 3] =    4.34;
  daa[ 9*20+ 4] =   62.73; daa[ 9*20+ 5] =    8.34; daa[ 9*20+ 6] =    3.31; daa[ 9*20+ 7] =    5.98;
  daa[ 9*20+ 8] =   12.26; daa[10*20+ 0] =   25.46; daa[10*20+ 1] =   15.58; daa[10*20+ 2] =   15.16;
  daa[10*20+ 3] =    1.90; daa[10*20+ 4] =   25.65; daa[10*20+ 5] =   39.70; daa[10*20+ 6] =    1.90;
  daa[10*20+ 7] =    2.41; daa[10*20+ 8] =   11.49; daa[10*20+ 9] =  329.09; daa[11*20+ 0] =    8.36;
  daa[11*20+ 1] =  141.40; daa[11*20+ 2] =  608.70; daa[11*20+ 3] =    2.31; daa[11*20+ 4] =    1.90;
  daa[11*20+ 5] =  465.58; daa[11*20+ 6] =  313.86; daa[11*20+ 7] =   22.73; daa[11*20+ 8] =  127.67;
  daa[11*20+ 9] =   19.57; daa[11*20+10] =   14.88; daa[12*20+ 0] =  141.88; daa[12*20+ 1] =    1.90;
  daa[12*20+ 2] =   65.41; daa[12*20+ 3] =    1.90; daa[12*20+ 4] =    6.18; daa[12*20+ 5] =   47.37;
  daa[12*20+ 6] =    1.90; daa[12*20+ 7] =    1.90; daa[12*20+ 8] =   11.97; daa[12*20+ 9] =  517.98;
  daa[12*20+10] =  537.53; daa[12*20+11] =   91.37; daa[13*20+ 0] =    6.37; daa[13*20+ 1] =    4.69;
  daa[13*20+ 2] =   15.20; daa[13*20+ 3] =    4.98; daa[13*20+ 4] =   70.80; daa[13*20+ 5] =   19.11;
  daa[13*20+ 6] =    2.67; daa[13*20+ 7] =    1.90; daa[13*20+ 8] =   48.16; daa[13*20+ 9] =   84.67;
  daa[13*20+10] =  216.06; daa[13*20+11] =    6.44; daa[13*20+12] =   90.82; daa[14*20+ 0] =   54.31;
  daa[14*20+ 1] =   23.64; daa[14*20+ 2] =   73.31; daa[14*20+ 3] =   13.43; daa[14*20+ 4] =   31.26;
  daa[14*20+ 5] =  137.29; daa[14*20+ 6] =   12.83; daa[14*20+ 7] =    1.90; daa[14*20+ 8] =   60.97;
  daa[14*20+ 9] =   20.63; daa[14*20+10] =   40.10; daa[14*20+11] =   50.10; daa[14*20+12] =   18.84;
  daa[14*20+13] =   17.31; daa[15*20+ 0] =  387.86; daa[15*20+ 1] =    6.04; daa[15*20+ 2] =  494.39;
  daa[15*20+ 3] =   69.02; daa[15*20+ 4] =  277.05; daa[15*20+ 5] =   54.11; daa[15*20+ 6] =   54.71;
  daa[15*20+ 7] =  125.93; daa[15*20+ 8] =   77.46; daa[15*20+ 9] =   47.70; daa[15*20+10] =   73.61;
  daa[15*20+11] =  105.79; daa[15*20+12] =  111.16; daa[15*20+13] =   64.29; daa[15*20+14] =  169.90;
  daa[16*20+ 0] =  480.72; daa[16*20+ 1] =    2.08; daa[16*20+ 2] =  238.46; daa[16*20+ 3] =   28.01;
  daa[16*20+ 4] =  179.97; daa[16*20+ 5] =   94.93; daa[16*20+ 6] =   14.82; daa[16*20+ 7] =   11.17;
  daa[16*20+ 8] =   44.78; daa[16*20+ 9] =  368.43; daa[16*20+10] =  126.40; daa[16*20+11] =  136.33;
  daa[16*20+12] =  528.17; daa[16*20+13] =   33.85; daa[16*20+14] =  128.22; daa[16*20+15] =  597.21;
  daa[17*20+ 0] =    1.90; daa[17*20+ 1] =   21.95; daa[17*20+ 2] =   10.68; daa[17*20+ 3] =   19.86;
  daa[17*20+ 4] =   33.60; daa[17*20+ 5] =    1.90; daa[17*20+ 6] =    1.90; daa[17*20+ 7] =   10.92;
  daa[17*20+ 8] =    7.08; daa[17*20+ 9] =    1.90; daa[17*20+10] =   32.44; daa[17*20+11] =   24.00;
  daa[17*20+12] =   21.71; daa[17*20+13] =    7.84; daa[17*20+14] =    4.21; daa[17*20+15] =   38.58;
  daa[17*20+16] =    9.99; daa[18*20+ 0] =    6.48; daa[18*20+ 1] =    1.90; daa[18*20+ 2] =  191.36;
  daa[18*20+ 3] =   21.21; daa[18*20+ 4] =  254.77; daa[18*20+ 5] =   38.82; daa[18*20+ 6] =   13.12;
  daa[18*20+ 7] =    3.21; daa[18*20+ 8] =  670.14; daa[18*20+ 9] =   25.01; daa[18*20+10] =   44.15;
  daa[18*20+11] =   51.17; daa[18*20+12] =   39.96; daa[18*20+13] =  465.58; daa[18*20+14] =   16.21;
  daa[18*20+15] =   64.92; daa[18*20+16] =   38.73; daa[18*20+17] =   26.25; daa[19*20+ 0] =  195.06;
  daa[19*20+ 1] =    7.64; daa[19*20+ 2] =    1.90; daa[19*20+ 3] =    1.90; daa[19*20+ 4] =    1.90;
  daa[19*20+ 5] =   19.00; daa[19*20+ 6] =   21.14; daa[19*20+ 7] =    2.53; daa[19*20+ 8] =    1.90;
  daa[19*20+ 9] = 1222.94; daa[19*20+10] =   91.67; daa[19*20+11] =    1.90; daa[19*20+12] =  387.54;
  daa[19*20+13] =    6.35; daa[19*20+14] =    8.23; daa[19*20+15] =    1.90; daa[19*20+16] =  204.54;
  daa[19*20+17] =    5.37; daa[19*20+18] =    1.90;
  
  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];

  pi[ 0] = 0.072000; pi[ 1] = 0.019000; pi[ 2] = 0.039000; pi[ 3] = 0.019000;
  pi[ 4] = 0.006000; pi[ 5] = 0.025000; pi[ 6] = 0.024000; pi[ 7] = 0.056000;
  pi[ 8] = 0.028000; pi[ 9] = 0.088000; pi[10] = 0.169000; pi[11] = 0.023000;
  pi[12] = 0.054000; pi[13] = 0.061000; pi[14] = 0.054000; pi[15] = 0.072000;
  pi[16] = 0.086000; pi[17] = 0.029000; pi[18] = 0.033000; pi[19] = 0.043000;

  return 1;
}


/*********************************************************/

model *Init_Model(allseq *data, option *input)
{
  /* Initialize substitution model */

  model *mod;
  int i,j;
  int ns;
  fit_double sum,aux;
  fit_double *di, *space;


  mod = input->mod;

  ns = mod->ns;

  mod->analytical = 1;

  For(i,6) mod->gtr_param[i] = 1.0;
  For(i,mod->n_catg)
    {
      mod->rr_mixturem[i] = 1.0;
      mod->r_proba[i]     = 1.0;
    }

  if(!input->mod->invar) For(i,data->crunch_len) data->invar[i] = 0;


  mod->update_eigen  = 0;

  space   = (fit_double *)mCalloc(2*ns,sizeof(fit_double));
  di      = (fit_double *)mCalloc(  ns,sizeof(fit_double));



  switch(mod->datatype)
    { 
    case NT :
      {

	switch(mod->model_applies_to)  
	  {
	  case NT :
	    {
	      Get_Base_Freqs(data,0,1);
    
	      For(i,ns) mod->pi[i] = data->b_frq[i];
	      
	      mod->lambda = 1.;
	      
	      switch(mod->subst_modelname)
		{
		case JC69 :
		  {
		    mod->kappa = 1.;
		    mod->pi[0] = mod->pi[1] = mod->pi[2] = mod->pi[3] = .25;
		    break;
		  }
		  
		case K80 :
		  {
		    mod->pi[0] = mod->pi[1] = mod->pi[2] = mod->pi[3] = .25;
		    break;
		  }
		  
		case F81 : 
		  {
		    mod->kappa = 1.;
		    break;
		  }
		  
		case  F84 :
		  {	      
		    aux = (mod->pi[0]+mod->pi[2]-mod->pi[1]-mod->pi[3])/(2.*mod->kappa);
		    mod->lambda = (mod->pi[1]+mod->pi[3] + aux)/(mod->pi[0]+mod->pi[2] -aux); 
		    break;
		  }
		  
		case GTR : 
		  {	  	  
		    mod->analytical   = 0;
		    mod->update_eigen = 1;
		    Update_Qmat_GTR(mod);
		    break;
		  }

		case HKY85 : break;

		default : 
		  {
		    char *s;
		    s = Return_Subst_Model_Name(mod);
		    printf("\n. Subst. model name : %s\n",s);
		    Free(s);
		    printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
		    Exit("");		    
		  }
		}

	      switch(mod->switch_modelname)
		{
		case NO_SWITCH : break;
		case SWITCH_S1 : break;
		case SWITCH_S2 : break;
		default : break;
		}
	      break;
	    }
	  case CODONS :
	    {
	      /* Init for codon model */
	     
	      Get_Codon_Freqs(data,mod);
	      
	      For(i,ns) mod->pi[i] = data->b_frq[i];
	      
	      mod->update_eigen  = 1;
	      mod->analytical    = 0;
	      
	      if(mod->switch_modelname == NO_SWITCH)
		{
		  int i;
		  
		  For(i,mod->n_catq)
		    {
		      mod->qmat_struct[0]->omega[i]             = mod->omega[i]; 
		      mod->qmat_struct[0]->omega_min[i]         = mod->omega_min[i]; 
		      mod->qmat_struct[0]->omega_max[i]         = mod->omega_max[i]; 
		      mod->qmat_struct[0]->omega_proba[i]       = 1.0;
		      mod->qmat_struct[0]->trans_omega_proba[i] = 1.0;
		      mod->qmat_struct[0]->qmat_proba[i]        = mod->omega_proba[i];
		      mod->qmat_struct[0]->trans_qmat_proba[i]  = mod->omega_proba[i];
		    }
		 
		  Init_T_Omega(mod);
		  
		  Init_T_Omega_Proba(mod);
		  
		  mod->qmat_struct[0]->curr_qmat_cat = 0;
		  
		  Update_Qmat_Codon(mod,mod->qmat_struct[0]);
		}
	      
	      else if((mod->switch_modelname == SWITCH_S1) ||
		      (mod->switch_modelname == SWITCH_S2))
		{
		  
		  /* if(mod->n_omega > 3) */
		  /*   { */
		  /*     printf("\n. Can't do more than 3 classes with switching models."); */
		  /*     printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
		  /*     Exit(""); */
		  /*   } */

		  mod->qmat_struct[0]->qmat_proba[0]       = 1.0;
		  mod->qmat_struct[0]->trans_qmat_proba[0] = 1.0;

		  For(i,mod->n_omega) 
		    {
		      mod->qmat_struct[0]->omega[i]             = mod->omega[i];
		      mod->qmat_struct[0]->omega_min[i]         = mod->omega_min[i];
		      mod->qmat_struct[0]->omega_max[i]         = mod->omega_max[i];
		      mod->qmat_struct[0]->omega_proba[i]       = mod->omega_proba[i];
		      mod->qmat_struct[0]->trans_omega_proba[i] = mod->omega_proba[i];
		    }

                  if(mod->io_delta_switch < 1.E-50) mod->io_delta_switch = 0.1;

		  For(i,mod->n_omega*(mod->n_omega-1)/2) 
		    mod->qmat_struct[0]->theta[i] = mod->io_delta_switch;
		  
		  /* mod->qmat_struct[0]->theta[0] = mod->io_delta_switch; */

		  Init_T_Omega(mod);
		  Init_T_Omega_Proba(mod);
		  		  		  
		  Update_Qmat_Codon(mod,mod->qmat_struct[0]);
		}
	      break;
	    }
	  default :  
	    {
	      printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
	      Exit("");
	    }
	  }
	break;
      }
    case AA : 
      { 
	fit_double *U, *V, *Root;
	
	mod->analytical = 0;
	
	U      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double));
	V      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double));
	Root   = (fit_double *)mCalloc(mod->ns,sizeof(fit_double));
	
	
	/* Init for amino-acids */
	
	Get_AA_Freqs(data); 
	
	switch(mod->subst_modelname)
	  {
	  case DAYHOFF : 
	    {
	      Init_Qmat_Dayhoff(mod->qmat_struct[0]->qmat,mod->pi);
	      break;
	    }
	  case JTT : 
	    {
	      Init_Qmat_JTT(mod->qmat_struct[0]->qmat, mod->pi);
	      break;
	    }
	  case MTREV : 
	    {
	      Init_Qmat_MtREV(mod->qmat_struct[0]->qmat, mod->pi);
	      break;
	    }
	  case WAG : 
	    {
	      Init_Qmat_WAG(mod->qmat_struct[0]->qmat, mod->pi);
	      break;
	    }
	  case DCMUT : 
	    {
	      Init_Qmat_DCMut(mod->qmat_struct[0]->qmat, mod->pi);
	      break;
	    }
	  default : break;
	  }


      For(i,ns) For(j,ns)
	mod->qmat_struct[0]->qmat[i*ns+j] *= mod->pi[j]/100.;
      
      
      mod->mr = .0;
      For (i,ns)
	{
	  sum=.0;
	  For(j, ns) sum += mod->qmat_struct[0]->qmat[i*ns+j];
	  mod->qmat_struct[0]->qmat[i*ns+i] = -sum;
	  mod->mr += mod->pi[i] * sum;
	}

      /* Scale instantaneous rate matrix */
      For (i,ns*ns) mod->qmat_struct[0]->qmat[i*ns+i] /= mod->mr;

      Update_Eigen(mod->ns,mod->qmat_struct[0]->qmat,U,V,Root,mod);
      
      For(i,mod->ns) For(j,mod->ns) 
	{
	  mod->qmat_struct[0]->u_mat[i*mod->ns+j] = U[i*mod->ns+j];
	  mod->qmat_struct[0]->v_mat[i*mod->ns+j] = V[i*mod->ns+j];
	}
      
      For(i,mod->ns)
	mod->qmat_struct[0]->expD_mr_vct[i] = (fit_double)exp(Root[i]);
      
  
      Free(U);
      Free(V);
      Free(Root);
      break;
      }
    default : 
      {
	printf("\n. %s %d\n",__FILE__,__LINE__);
	Exit("\n. Err. \n");
      }
    }
  free(di);free(space);
  return mod;
}

/*********************************************************/

int PMat_Numeric(fit_double l, int ns, fit_double **Pij, qmat *qmat_struct, arbre *tree)
{
  /* Compute change probabilities using eigen vectors and eigen values */

  int n;
  int i, j, k;
  fit_double *U,*V,*R;
  fit_double *expt;
  fit_double *uexpt;
  int dim1,dim2,neg;
  fit_double sum;

  dim2 = ns*ns;
  dim1 = ns;


  n = ns;

  expt  = (fit_double *)mCalloc(n,  sizeof(fit_double));
  uexpt = (fit_double *)mCalloc(n*n,sizeof(fit_double));


  U = qmat_struct->u_mat       + dim2*qmat_struct->curr_qmat_cat;
  V = qmat_struct->v_mat       + dim2*qmat_struct->curr_qmat_cat;
  R = qmat_struct->expD_mr_vct + dim1*qmat_struct->curr_qmat_cat;

  /* Compute pow((fit_double)exp(D/mr),l) into mat_eDmrl */
  For(k,n) expt[k] = (fit_double)pow(R[k],l);

/*   For(i,n) For(k,n) printf("%f %f\n",qmat_struct->curr_qmat_cat,U[i*n+k]); */
  /* Multiply Vr*pow((fit_double)exp(D/mr),l)*Vi into Pij */
  For(i,n) For(k,n) uexpt[i*n+k] = U[i*n+k] * expt[k];
  
  neg = 0;
  For(i,n)
    {      
      sum = .0;
      For(j,n)
	{
	  Pij[i][j] = .0;
	  
	  For(k,n)  Pij[i][j] += uexpt[i*n+k] * V[k*n+j];
	  
	  /* if(Pij[i][j] < MDBL_MIN) */
	  /*   { */
	  /*     Pij[i][j] = 1.E-30; */
	  /*   } */
	  if(Pij[i][j] < 0.0)
	    {
	      if(Pij[i][j] < 0.0-1.E-5) neg = 1;
	      /* Pij[i][j] = 1.E-60; */
	      Pij[i][j] = 0.0;
	    }
	  sum += Pij[i][j];
	  
	}
      if((sum < 1.0-1.E-4) || (sum > 1.0+1.E-4))
	{
          /* FILE *fp_err; */
          /* fp_err = fopen("/tmp/fitmodel.err","w"); */
          /* sum = .0; */
          
          /* For(j,n)  */
          /*   { */
          /*     sum += Pij[i][j]; */
          /*     fprintf(fp_err,"[Pij %3d - %3d]  %G %G\n",i,j,Pij[i][j],sum); */
          /*   } */

          /* For(i,n) For(k,n) fprintf(fp_err,"[UV %3d - %3d]  %f %f %G\n",i,k,U[i*n+k],V[i*n+k],uexpt[i*n+k]); */
          /* For(i,n) fprintf(fp_err,"[L  %3d]  %G %G %G\n",i,expt[i],R[i],l); */
          /* fclose(fp_err); */

          printf("\n= WARNING: numerical precision issue encountered when computing Pij matrix.");
	  printf("\n= Sum of probas = %f\n",sum);
          Print_Param(stderr,tree);
	  /* Exit("\n"); */

          sum = 0.0;
          For(j,n) sum += Pij[i][j];
          For(j,n) Pij[i][j] /= sum;
	}
    }
  Free(expt); Free(uexpt);
  if(neg) 
    {
      printf("neg encountered  l=%f qmatcat = %d\n",
	     l,
	     qmat_struct->curr_qmat_cat);
      fflush(NULL);
      return 0;
    }
  return 1;
}

/* /\*********************************************************\/ */

void Update_Qmat_GTR(model *mod)
{
  /* Fill the instantaneous rate matrix for the GTR model */

  int i,j;
  fit_double *U, *V, *Root;


  U      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double));
  V      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double));
  Root   = (fit_double *)mCalloc(mod->ns,sizeof(fit_double));

  if((mod->gtr_param[5] < 0.9999) || (mod->gtr_param[5] > 1.0001))
    Exit("\n. Err. in Update_Qmat_GTR\n");
  
  mod->qmat_struct[0]->qmat[0*4+1] = mod->gtr_param[0]*mod->pi[1];
  mod->qmat_struct[0]->qmat[0*4+2] = mod->gtr_param[1]*mod->pi[2];
  mod->qmat_struct[0]->qmat[0*4+3] = mod->gtr_param[2]*mod->pi[3];
  
  mod->qmat_struct[0]->qmat[1*4+0] = mod->gtr_param[0]*mod->pi[0];
  mod->qmat_struct[0]->qmat[1*4+2] = mod->gtr_param[3]*mod->pi[2];
  mod->qmat_struct[0]->qmat[1*4+3] = mod->gtr_param[4]*mod->pi[3];
  
  mod->qmat_struct[0]->qmat[2*4+0] = mod->gtr_param[1]*mod->pi[0];
  mod->qmat_struct[0]->qmat[2*4+1] = mod->gtr_param[3]*mod->pi[1];
  mod->qmat_struct[0]->qmat[2*4+3] = mod->gtr_param[5]*mod->pi[3];
  
  mod->qmat_struct[0]->qmat[3*4+0] = mod->gtr_param[2]*mod->pi[0];
  mod->qmat_struct[0]->qmat[3*4+1] = mod->gtr_param[4]*mod->pi[1];
  mod->qmat_struct[0]->qmat[3*4+2] = mod->gtr_param[5]*mod->pi[2];

  mod->qmat_struct[0]->qmat[0*4+0] = -(mod->gtr_param[0]*mod->pi[1]+mod->gtr_param[1]*mod->pi[2]+mod->gtr_param[2]*mod->pi[3]);
  mod->qmat_struct[0]->qmat[1*4+1] = -(mod->gtr_param[0]*mod->pi[0]+mod->gtr_param[3]*mod->pi[2]+mod->gtr_param[4]*mod->pi[3]);
  mod->qmat_struct[0]->qmat[2*4+2] = -(mod->gtr_param[1]*mod->pi[0]+mod->gtr_param[3]*mod->pi[1]+mod->gtr_param[5]*mod->pi[3]);
  mod->qmat_struct[0]->qmat[3*4+3] = -(mod->gtr_param[2]*mod->pi[0]+mod->gtr_param[4]*mod->pi[1]+mod->gtr_param[5]*mod->pi[2]);

  mod->mr = .0;
  For(i,mod->ns) mod->mr -= mod->pi[i] * (mod->qmat_struct[0]->qmat[i*mod->ns+i]);
  For(i,mod->ns*mod->ns) mod->qmat_struct[0]->qmat[i] /= mod->mr;


  Update_Eigen(mod->ns,mod->qmat_struct[0]->qmat,U,V,Root,mod);
  
  For(i,mod->ns) For(j,mod->ns)
    {
      mod->qmat_struct[0]->u_mat[i*mod->ns+j] = U[i*mod->ns+j];
      mod->qmat_struct[0]->v_mat[i*mod->ns+j] = V[i*mod->ns+j];
    }
  
  For(i,mod->ns)
    mod->qmat_struct[0]->expD_mr_vct[i] = (fit_double)exp(Root[i]);

  
  Free(U);
  Free(V);
  Free(Root);

}

/* /\* /\\*********************************************************\\/ *\/ */

/* void Update_Qmat_GTR(model *mod) */
/* { */
/*   /\* Fill the instantaneous rate matrix for the GTR model *\/ */

/*   int i,j; */
/*   fit_double *U, *V, *Root; */
/*   qmat *qmat_struct; */
/*   int class_i, class_j; */

/*   U      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double)); */
/*   V      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double)); */
/*   Root   = (fit_double *)mCalloc(mod->ns,        sizeof(fit_double)); */

/*   if((mod->gtr_param[5] < 0.9999) || (mod->gtr_param[5] > 1.0001)) */
/*     Exit("\n. Err. in Update_Qmat_GTR\n"); */
  
/*   qmat_struct = mod->qmat_struct[0]; */

/*   For(i,mod->ns*mod->ns) qmat_struct->qmat[i] = .0; */

/*   For(i,mod->ns) */
/*     { */
/*       For(j,mod->ns) */
/* 	{ */
/* 	  if(i != j) /\* We'll fill the diagonal later *\/ */
/* 	    { */
/* 	      if((i/mod->ns_nt) == (j/mod->ns_nt)) /\* substitutions btw nucleotides *\/ */
/* 		{	       */
/* 		  qmat_struct->qmat[i*mod->ns_nt+j] =  */
/* /\* 		    Return_Relative_Rate_GTR(i%mod->ns_nt,j%mod->ns_nt,mod) * *\/ */
/* 		    mod->gtr_param[MMIN(i%mod->ns_nt,j%mod->ns_nt)*mod->ns_nt+MMAX(j%mod->ns_nt,i%mod->ns_nt) -  */
/* 				  (MMIN(i%mod->ns_nt,j%mod->ns_nt)+1+(int)pow(MMIN(i%mod->ns_nt,j%mod->ns_nt)+1,2))/2] *  */
/* 		    mod->rsubst_rate[(int)i/mod->ns_nt]; */
/* 		} */
/* 	      else /\* switches btw classes *\/ */
/* 		{ */
/* 		  class_i = (int)i/mod->ns_nt; */
/* 		  class_j = (int)j/mod->ns_nt; */

/* 		  if((i%mod->ns_nt) == (j%mod->ns_nt)) */
/* 		    qmat_struct->qmat[i*mod->ns_nt+j] =  */
/* 		      qmat_struct->theta[MMIN(class_i,class_j)*mod->ns_nt+MMAX(class_j,class_i) -  */
/* 					(MMIN(class_i,class_j)+1+(int)pow(MMIN(class_i,class_j)+1,2))/2]; */
/* 		} */
/* 	    } */
/* 	} */
/*     } */





/*   Update_Eigen(mod->ns,qmat_struct->qmat,U,V,Root,mod); */
  
/*   For(i,mod->ns) For(j,mod->ns)  */
/*     { */
/*       qmat_struct->u_mat[i*mod->ns+j] = U[i*mod->ns+j]; */
/*       qmat_struct->v_mat[i*mod->ns+j] = V[i*mod->ns+j]; */
/*     } */
  
/*   For(i,mod->ns) */
/*     qmat_struct->expD_mr_vct[i] = (fit_double)exp(Root[i]); */

  
/*   Free(U); */
/*   Free(V); */
/*   Free(Root); */

/* } */

/*********************************************************/

void Set_Model_Parameters(arbre *tree)
{
  /* Set the parameters of the substitution model before the eigenvalue decomposition */

    fit_double sum;
    int i;

    DiscreteGamma (tree->mod->r_proba, tree->mod->rr_mixturem, tree->mod->alpha,
                   tree->mod->alpha,tree->mod->n_catg,0);
    
    if(tree->mod->model_applies_to == CODONS)
      { 


	if(tree->mod->switch_modelname == NO_SWITCH) 
	  {

            For(i,tree->mod->n_catq)
              {
                if(tree->mod->qmat_struct[0]->omega[i] < tree->mod->qmat_struct[0]->omega_min[i])
                  tree->mod->qmat_struct[0]->omega[i] = tree->mod->qmat_struct[0]->omega_min[i];
                if(tree->mod->qmat_struct[0]->omega[i] > tree->mod->qmat_struct[0]->omega_max[i])
                  tree->mod->qmat_struct[0]->omega[i] = tree->mod->qmat_struct[0]->omega_max[i];
              }

	    Untransform_Probs(tree->mod->qmat_struct[0]->trans_qmat_proba,
			      tree->mod->qmat_struct[0]->qmat_proba,
			      tree->mod->n_catq);
	    
	  }
	else if((tree->mod->switch_modelname == SWITCH_S1) || (tree->mod->switch_modelname == SWITCH_S2))
	  {
	    
            For(i,tree->mod->n_omega)
              {
                if(tree->mod->qmat_struct[0]->omega[i] < tree->mod->qmat_struct[0]->omega_min[i])
                  tree->mod->qmat_struct[0]->omega[i] = tree->mod->qmat_struct[0]->omega_min[i];
                if(tree->mod->qmat_struct[0]->omega[i] > tree->mod->qmat_struct[0]->omega_max[i])
                  tree->mod->qmat_struct[0]->omega[i] = tree->mod->qmat_struct[0]->omega_max[i];
              }

	    Untransform_Probs(tree->mod->qmat_struct[0]->trans_omega_proba,
			      tree->mod->qmat_struct[0]->omega_proba,
			      tree->mod->n_omega);
	    
	    For(i,tree->mod->n_omega)
		tree->mod->qmat_struct[0]->omega[i] = (fit_double)fabs(tree->mod->qmat_struct[0]->omega[i]);

	    For(i,tree->mod->n_omega*(tree->mod->n_omega-1)/2)
	      tree->mod->qmat_struct[0]->theta[i] = (fit_double)fabs(tree->mod->qmat_struct[0]->theta[i]);

	  }
	else
	  {
	    printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
	    Exit("");		    
	  }
      }
    
	
    if(tree->mod->model_applies_to == CODONS)
      {
	if((tree->mod->switch_modelname == SWITCH_S1) || (tree->mod->switch_modelname == SWITCH_S2))
	  {
	    if(tree->mod->switch_modelname == SWITCH_S1) /* S1 */
	      {
		for(i=1;i<tree->mod->n_omega*(tree->mod->n_omega-1)/2;i++)
		  tree->mod->qmat_struct[0]->theta[i] = tree->mod->qmat_struct[0]->theta[0];
	      }
	    For(i,tree->mod->n_omega*(tree->mod->n_omega-1)/2)
	      {
		if(tree->mod->qmat_struct[0]->theta[i] > 500.0) tree->mod->qmat_struct[0]->theta[i] = 500.0;
		if(tree->mod->qmat_struct[0]->theta[i] < 1.E-3) tree->mod->qmat_struct[0]->theta[i] = 1.E-3;
	      }
	  }
	
	if(tree->mod->switch_modelname == NO_SWITCH)
	  {
	    For(i,tree->mod->n_omega) tree->mod->qmat_struct[0]->omega_proba[i] = 1.0;
	  }
      }
    
    else if((tree->mod->model_number < 10) && (tree->mod->s_opt->opt_bfreq))
      {
	sum = .0;
	For(i,4) sum += tree->mod->pi[i];
	For(i,4) 
	  {
	    tree->mod->pi[i] /= sum;
	    /* 	  printf("pi[%d]->%f\n",i+1,tree->mod->pi[i]); */
	  }
      }
    
    if(tree->mod->subst_modelname == GTR)
      {
	For(i,6) tree->mod->gtr_param[i] = (fit_double)fabs(tree->mod->gtr_param[i]);
 	For(i,6) 
	  tree->mod->gtr_param[i] /= tree->mod->gtr_param[5];
      }
    
    if(tree->mod->update_eigen) 
      {
	if(tree->mod->model_applies_to == NT)          Update_Qmat_GTR(tree->mod);
	else if(tree->mod->model_applies_to == CODONS) Update_Qmat_Codon(tree->mod,tree->mod->qmat_struct[0]);
      }
    
    
    /* Only used with serial sample data */
/*     if(tree->mod->tpos_ols) */
/*       Tpos_OLS(tree,tree->root,tree->root->v[0]); */
    
/*     if(tree->mod->update_bl_using_tpos) */
/*         { */
/*             tree->root->tpos = .0; */
/*             Update_Bl_Using_Tpos(tree,tree->root,tree->root->v[0]); */
/*         } */
    
/*     if(tree->mod->update_bl_using_tpos || tree->mod->tpos_ols) */
/*         Update_Tpos_Min_Post(tree->root,tree->root->v[0]); */
}

/*********************************************************/

/* void Update_Qmat_Codon(model *mod, qmat *qmat_struct) */
/* { */
/*   /\* Fill the instantaneous rate matrix (codon-based model */
/*      of substitutions) */
/*   *\/ */

/*   int i,j,k; */
/*   int aa1,aa2; */
/*   fit_double sum; */
/*   int num_codon1, num_codon2; */
/*   int dim2,dim3; */
/*   int n_catq; */
/*   fit_double *w; */
/*   fit_double *U, *V, *Root; */
/*   fit_double div; */
/*   int codon_64_i, codon_64_j; */
/*   int cod_i, cod_j; */

/*   div = 0.0; */

/*   U      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double)); */
/*   V      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double)); */
/*   Root   = (fit_double *)mCalloc(mod->ns,sizeof(fit_double)); */
 
/*   w = qmat_struct->omega; */

/*   n_catq   = mod->n_catq; */
/*   dim2     = mod->ns; */
/*   dim3     = mod->ns*mod->ns; */

/*   num_codon1 = num_codon2 = -1; */


/*   ////////////////////////// */

/* /\*   qmat_struct->qmat_proba[0]         = 0.66; *\/ */
/* /\*   qmat_struct->qmat_proba[1]         = 0.30; *\/ */
/* /\*   qmat_struct->qmat_proba[2]         = 0.04; *\/ */
/* /\*   (*(qmat_struct->t_omega[0]))       = 0.00; *\/ */
/* /\*   (*(qmat_struct->t_omega[3]))       = 1.00; *\/ */
/* /\*   (*(qmat_struct->t_omega[6]))       = 9.72; *\/ */

/* /\*   mod->kappa                         = 10.401; *\/ */


/* /\*   (*(qmat_struct->t_omega_proba[0])) = 0.66; *\/ */
/* /\*   (*(qmat_struct->t_omega_proba[1])) = 0.30; *\/ */
/* /\*   (*(qmat_struct->t_omega_proba[2])) = 0.04; *\/ */
/* /\*   (*(qmat_struct->t_omega[0]))       = 0.00; *\/ */
/* /\*   (*(qmat_struct->t_omega[1]))       = 1.00; *\/ */
/* /\*   (*(qmat_struct->t_omega[2]))       = 9.72; *\/ */
/* /\*   qmat_struct->theta[0]              = 0.0001; *\/ */
/* /\*   qmat_struct->theta[1]              = 0.0001; *\/ */
/* /\*   qmat_struct->theta[2]              = 0.0001; *\/ */

  
/*   ////////////////////////// */

/*   mod->delta_switch = .0; */
/*   if(mod->switch_modelname != NO_SWITCH) */
/*     { */
/*       For(i,mod->n_omega) */
/* 	{ */
/* 	  For(j,mod->n_omega) */
/* 	    { */
/* 	      if(i != j) */
/* 		{ */
/* 		  mod->delta_switch += */
/* 		    (*(qmat_struct->t_omega_proba[i])) * */
/* 		    ((*(qmat_struct->t_omega_proba[j])) * */
/* 		     qmat_struct->theta[MMIN(i,j) * mod->n_omega + MMAX(i,j) - */
/* 				       (MMIN(i,j)+1+(int)pow(MMIN(i,j)+1,2))/2]); */
/* 		} */
/* 	    } */
/* 	} */
/*     } */
 
/*   For(i,mod->ns) */
/*     { */
/*       codon_64_i = i%mod->ns_codon; */
/*       num_codon1 = mod->c_code->sense_c[codon_64_i]; */

/*       for(j=i;j<mod->ns;j++) */
/* 	{ */
/* 	  codon_64_j = j%mod->ns_codon; */
/* 	  num_codon2 = mod->c_code->sense_c[codon_64_j]; */

/* 	  if((i/mod->ns_codon) != (j/mod->ns_codon)) */
/* 	     { */
/* 	       if((i%mod->ns_codon) == (j%mod->ns_codon)) */
/* 		 { */
/* 		   For(k,n_catq) */
/* 		     { */
/* 		       cod_i = (int)i/mod->ns_codon; */
/* 		       cod_j = (int)j/mod->ns_codon; */

/* 		       qmat_struct->qmat[k*dim3+i*mod->ns+j] = */
/* 			 qmat_struct->theta[MMIN(cod_i,cod_j) * mod->n_omega + MMAX(cod_i,cod_j) - */
/* 					   (MMIN(cod_i,cod_j)+1+(int)pow(MMIN(cod_i,cod_j)+1,2))/2]; */


/* 		       qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j]; */


/* 		       if(qmat_struct->qmat[k*dim3+i*mod->ns+j] < 0.0) */
/* 			 { */
/* 			   printf("%f %f\n", */
/* 				  qmat_struct->theta[MMIN(cod_i,cod_j) * mod->n_omega + MMAX(cod_i,cod_j) - */
/* 						    (MMIN(cod_i,cod_j)+1+(int)pow(MMIN(cod_i,cod_j)+1,2))/2], */
/* 				  qmat_struct->qmat[k*dim3+i*mod->ns+j]); */
/* 			   Exit("\n. Q(i!=j) < 0.0 !?!\n"); */
/* 			 } */
/* 		     } */
/* 		 } */
/* 	       else */
/* 		 { */
/* 		   For(k,n_catq) */
/* 		     { */
/* 		       qmat_struct->qmat[k*dim3+i*mod->ns+j] = 0.0; */
/* 		       qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j]; */
/* 		     } */
/* 		 } */
/* 	     } */
/* 	  else */
/* 	    { */
/* 	      if(mod->c_code->n_diff_b_2_codons[codon_64_i*mod->ns_codon+codon_64_j] == 1) */
/* 		{ */
/* 		  aa1 = Get_AA(mod->c_code,num_codon1); */
/* 		  aa2 = Get_AA(mod->c_code,num_codon2); */
		  
/* 		  For(k,n_catq) */
/* 		    { */
/* 		      qmat_struct->qmat[k*dim3+i*mod->ns+j]= */
/* 			(aa1 == aa2)? */
/* 			( */
/* 			 (mod->c_code->tstvtable[codon_64_i*mod->ns_codon+codon_64_j]==1)? */
/* 			 (mod->kappa): */
/* 			 (1.0) */
/* 			 ): */
/* 			( */
/* 			 (mod->c_code->tstvtable[codon_64_i*mod->ns_codon+codon_64_j]==1)? */
/* 			 (*(qmat_struct->t_omega[k*mod->n_omega+ */
/* 						 (int)i/mod->ns_codon])*mod->kappa): */
/* 			 (*(qmat_struct->t_omega[k*mod->n_omega+ */
/*  						 (int)i/mod->ns_codon])) */
/* 			 ); */


/* 		      qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j]; */

/* 		      if(qmat_struct->qmat[k*dim3+i*mod->ns+j] < 0.0) */
/* 			{ */
/* 			  printf("%f %f %f\n", */
/* 				 mod->kappa, */
/* 				 *(qmat_struct->t_omega[k*mod->n_omega+ */
/* 							(int)i/mod->ns_codon]), */
/* 				 qmat_struct->qmat[k*dim3+i*mod->ns+j]); */
/* 			  Exit("\n. Q(i!=j) < 0.0 !?!\n"); */
/* 			} */


/* 		    } */
/* 		} */
/* 	      else */
/* 		For(k,n_catq) */
/* 		  { */
/* 		    qmat_struct->qmat[k*dim3+i*dim2+j] = 0.0; */
/* 		    qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j]; */
/* 		  } */
/* 	    } */
/* 	} */
/*     } */



/*   /\* Build equilibrium frequencies into the rate matrix *\/ */
/*   For(k,n_catq) */
/*     { */
/*       For(i,mod->ns) */
/* 	{ */
/* 	  qmat_struct->qmat[k*dim3+i*mod->ns+i] = 0.0; */
               
/* 	  For(j, mod->ns) */
/* 	    { */
/* 	      if((i/mod->ns_codon) != (j/mod->ns_codon)) */
/* 		{ */
/* 		  qmat_struct->qmat[k*dim3+i*dim2+j] *= */
/* 		    (*(qmat_struct->t_omega_proba[k*mod->n_omega + (int)j/mod->ns_codon])); */
/* 		} */
/* 	      else */
/* 		{ */
/* 		  qmat_struct->qmat[k*dim3+i*dim2+j] *= */
/* 		    mod->pi[j%mod->ns_codon]; */
/* 		} */
/* 	    } */
/* 	} */
/*     } */

/*   /\* Scaling factor so that the expected number of codon substitutions per time unit = 1 *\/ */
/*   mod->mr = .0; */
/*   For(k,n_catq) */
/*     { */
/*       For(i,mod->ns) */
/* 	{ */
/* 	  For(j,mod->ns) */
/* 	    { */
/* 	      if((int)(i/mod->ns_codon) == (int)(j/mod->ns_codon)) */
/* 		{ */
/* 		  mod->mr += */
/* 		    mod->pi[i%mod->ns_codon] * */
/* 		    (*(qmat_struct->t_omega_proba[(int)i/mod->ns_codon])) * */
/* 		    qmat_struct->qmat_proba[k] * */
/* 		    qmat_struct->qmat[k*dim3+i*dim2+j]; */
/* 		} */
/* 	    } */
/* 	} */
/*     } */

/* /\*     printf("SCALER = %f\n",mod->mr); *\/ */


/*   /\* Scale matrix *\/ */
/*   For(k,n_catq) */
/*     { */
/*       For(i,mod->ns) */
/* 	{ */
/* 	  For(j,mod->ns) */
/* 	    { */
/* 	      if((int)(i/mod->ns_codon) == (int)(j/mod->ns_codon)) */
/* 		{ */
/* 		  qmat_struct->qmat[k*dim3+i*dim2+j] /= mod->mr; */
/* 		} */
/* 	    } */
/* 	} */
/*     } */



/*   /\* Diagonal  terms *\/ */
/*   For(k,n_catq) */
/*     { */
/*       For(i,mod->ns) */
/* 	{ */
/* 	  sum = .0; */
/* 	  For(j,mod->ns) */
/* 	    sum += qmat_struct->qmat[k*dim3+i*dim2+j]; */
/* 	  qmat_struct->qmat[k*dim3+i*dim2+i] = -sum; */
/* 	} */
/*     } */

/*   For(k,n_catq) */
/*     { */

/*       Update_Eigen(mod->ns,qmat_struct->qmat+k*dim3,U,V,Root,mod); */

/*       For(i,mod->ns) For(j,mod->ns) */
/* 	{ */
/* 	  qmat_struct->u_mat[k*dim3+i*dim2+j]  = U[i*dim2+j]; */
/* 	  qmat_struct->v_mat[k*dim3+i*dim2+j]  = V[i*dim2+j]; */
/* 	  if(V[i*dim2+j] > 1.E+20) */
/* 	    Exit("\n. Numerical overflow during the calculation of eigen vectors\n"); */
/* 	} */

/*       For(i,mod->ns) */
/* 	{ */
/* 	  qmat_struct->expD_mr_vct[k*dim2+i] = (fit_double)exp(Root[i]); */
/* 	} */
/*     } */

/*   Free(U); */
/*   Free(V); */
/*   Free(Root); */
/* } */

/*********************************************************/

void Update_Qmat_Codon(model *mod, qmat *qmat_struct)
{
  /* Fill the instantaneous rate matrix (codon-based model
     of substitutions)
  */
  
  int i,j,k;
  int aa1,aa2;
  fit_double sum;
  int num_codon1, num_codon2;
  int dim2,dim3;
  int n_catq;
  fit_double *U, *V, *Root;
  int codon_64_i, codon_64_j;
  int cod_i, cod_j;
  
  
  U      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double));
  V      = (fit_double *)mCalloc(mod->ns*mod->ns,sizeof(fit_double));
  Root   = (fit_double *)mCalloc(mod->ns,sizeof(fit_double));
  
  
  n_catq   = mod->n_catq;
  dim2     = mod->ns;
  dim3     = mod->ns*mod->ns;
  
  num_codon1 = num_codon2 = -1;
  

  //////////////////////////

/*   qmat_struct->qmat_proba[0]         = 0.66; */
/*   qmat_struct->qmat_proba[1]         = 0.30; */
/*   qmat_struct->qmat_proba[2]         = 0.04; */
/*   (*(qmat_struct->t_omega[0]))       = 0.00; */
/*   (*(qmat_struct->t_omega[3]))       = 1.00; */
/*   (*(qmat_struct->t_omega[6]))       = 9.72; */

/*   mod->kappa                         = 10.401; */


/*   (*(qmat_struct->t_omega_proba[0])) = 0.66; */
/*   (*(qmat_struct->t_omega_proba[1])) = 0.30; */
/*   (*(qmat_struct->t_omega_proba[2])) = 0.04; */
/*   (*(qmat_struct->t_omega[0]))       = 0.00; */
/*   (*(qmat_struct->t_omega[1]))       = 1.00; */
/*   (*(qmat_struct->t_omega[2]))       = 9.72; */

/*   qmat_struct->theta[0]              = 0.0001; */
/*   qmat_struct->theta[1]              = 0.0001; */
/*   qmat_struct->theta[2]              = 0.0001; */

  
  //////////////////////////

  /* This loop only concerns codon substitutions */
  For(i,mod->ns)
    {
      codon_64_i = i%mod->ns_codon;
      num_codon1 = mod->c_code->sense_c[codon_64_i];

      for(j=i;j<mod->ns;j++)
	{
	  codon_64_j = j%mod->ns_codon;
	  num_codon2 = mod->c_code->sense_c[codon_64_j];

	  if((int)(i/mod->ns_codon) == (int)(j/mod->ns_codon))
	    {
	      if(mod->c_code->n_diff_b_2_codons[codon_64_i*mod->ns_codon+codon_64_j] == 1)
		{
		  aa1 = Get_AA(mod->c_code,num_codon1);
		  aa2 = Get_AA(mod->c_code,num_codon2);
		  
		  For(k,n_catq)
		    {
		      qmat_struct->qmat[k*dim3+i*mod->ns+j]=
			(aa1 == aa2)?
			(
			 (mod->c_code->tstvtable[codon_64_i*mod->ns_codon+codon_64_j]==1)?
			 (mod->kappa):
			 (1.0)
			 ):
			(
			 (mod->c_code->tstvtable[codon_64_i*mod->ns_codon+codon_64_j]==1)?
			 (*(qmat_struct->t_omega[k*mod->n_omega+
						 (int)i/mod->ns_codon])*mod->kappa):
			 (*(qmat_struct->t_omega[k*mod->n_omega+
 						 (int)i/mod->ns_codon]))
			 );


		      qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j];
		      
		      if(qmat_struct->qmat[k*dim3+i*mod->ns+j] < 0.0)
			{
			  printf("%f %f %f\n",
				 mod->kappa,
				 *(qmat_struct->t_omega[k*mod->n_omega+
							(int)i/mod->ns_codon]),
				 qmat_struct->qmat[k*dim3+i*mod->ns+j]);
			  Exit("\n. Q(i!=j) < 0.0 !?!\n");
			}
		    }
		}
	      else
		For(k,n_catq)
		  {
		    qmat_struct->qmat[k*dim3+i*dim2+j]    = 0.0;
		    qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j];
		  }
	    }
	}
    }

  /* Build equilibrium frequencies into the rate matrix */
  For(k,n_catq)
    {
      For(i,mod->ns)
	{
	  qmat_struct->qmat[k*dim3+i*mod->ns+i] = 0.0;
	  
	  For(j, mod->ns)
	    {
	      if((int)(i/mod->ns_codon) == (int)(j/mod->ns_codon))
		{
		  qmat_struct->qmat[k*dim3+i*dim2+j] *= mod->pi[j%mod->ns_codon];
		}
	    }
	}
    }

  /* Scaling factor = expected number of codon substitutions per time unit */
  mod->mr = .0;
  For(k,n_catq)
    {
      For(i,mod->ns)
	{
	  For(j,mod->ns)
	    {
	      if((int)(i/mod->ns_codon) == (int)(j/mod->ns_codon))
		{
		  mod->mr +=
		    mod->pi[i%mod->ns_codon] *
		    (*(qmat_struct->t_omega_proba[(int)i/mod->ns_codon])) *
		    qmat_struct->qmat_proba[k] *
		    qmat_struct->qmat[k*dim3+i*dim2+j];
		}
	    }
	}
    }

  if(mod->switch_modelname != NO_SWITCH)
    {
      mod->delta_switch = .0;

      For(i,mod->n_omega)
	{
	  For(j,mod->n_omega)
	    {
	      if(i != j)
		{
/* 		  mod->delta_switch += */
/* 		    (*(qmat_struct->t_omega_proba[i])) * */
/* 		    ((*(qmat_struct->t_omega_proba[j])) * */
/* 		     qmat_struct->theta[MMIN(i,j) * mod->n_omega + MMAX(i,j) - */
/* 				       (MMIN(i,j)+1+(int)pow(MMIN(i,j)+1,2))/2])/mod->mr; */
		  mod->delta_switch +=
		    (*(qmat_struct->t_omega_proba[i])) *
		    (*(qmat_struct->t_omega_proba[j])) *
		    qmat_struct->theta[MMIN(i,j) * mod->n_omega + MMAX(i,j) -
				      (MMIN(i,j)+1+(int)pow(MMIN(i,j)+1,2))/2];
		}
	    }
	}
      
      For(i,mod->ns)
	{
	  codon_64_i = i%mod->ns_codon;
	  num_codon1 = mod->c_code->sense_c[codon_64_i];
	  
	  for(j=i;j<mod->ns;j++)
	    {
	      codon_64_j = j%mod->ns_codon;
	      num_codon2 = mod->c_code->sense_c[codon_64_j];
	      
	      if((int)(i/mod->ns_codon) != (int)(j/mod->ns_codon))
		{
		  if((int)(i%mod->ns_codon) == (int)(j%mod->ns_codon))
		    {
		      For(k,n_catq)
			{
			  cod_i = (int)i/mod->ns_codon;
			  cod_j = (int)j/mod->ns_codon;
			  
			  qmat_struct->qmat[k*dim3+i*mod->ns+j] =
			    qmat_struct->theta[MMIN(cod_i,cod_j) * mod->n_omega + MMAX(cod_i,cod_j) -
					      (MMIN(cod_i,cod_j)+1+(int)pow(MMIN(cod_i,cod_j)+1,2))/2];
			  
			  qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j];
			  			  
			  if(qmat_struct->qmat[k*dim3+i*mod->ns+j] < 0.0)
			    {
			      printf("%f %f\n",
				     qmat_struct->theta[MMIN(cod_i,cod_j) * mod->n_omega +
							MMAX(cod_i,cod_j) -
							(MMIN(cod_i,cod_j)+1+
							 (int)pow(MMIN(cod_i,cod_j)+1,2))/2],
				     qmat_struct->qmat[k*dim3+i*mod->ns+j]);
			      Exit("\n. Q(i!=j) < 0.0 !?!\n");
			    }
			}
		    }
		  else
		    {
		      For(k,n_catq)
			{
			  qmat_struct->qmat[k*dim3+i*mod->ns+j] = 0.0;
			  qmat_struct->qmat[k*dim3+j*mod->ns+i] = qmat_struct->qmat[k*dim3+i*mod->ns+j];
			}
		    }
		}
	    }
	}


      /* Build equilibrium frequencies into the rate matrix */
      For(i,mod->ns)
	{
	  For(j,mod->ns)
	    {
	      if((int)(i/mod->ns_codon) != (int)(j/mod->ns_codon))
		{
		  qmat_struct->qmat[i*mod->ns+j] *=
		    (*(qmat_struct->t_omega_proba[(int)j/mod->ns_codon]));
		}
	    }
	}       
    }


  /* Scale the matrix */
  For(k,n_catq)
    {
      For(i,mod->ns)
	{
	  For(j,mod->ns)
	    {
	      if((int)(i/mod->ns_codon) == (int)(j/mod->ns_codon))
		qmat_struct->qmat[k*dim3+i*dim2+j] /= mod->mr;
	    }
	}
    }


  /* Diagonal  terms */
  For(k,n_catq)
    {
      For(i,mod->ns)
	{
	  sum = .0;
	  For(j,mod->ns)
	    sum += qmat_struct->qmat[k*dim3+i*dim2+j];
	  qmat_struct->qmat[k*dim3+i*dim2+i] = -sum;
	}
    }



  For(k,n_catq)
    {

      Update_Eigen(mod->ns,qmat_struct->qmat+k*dim3,U,V,Root,mod);

      For(i,mod->ns) For(j,mod->ns)
	{
	  qmat_struct->u_mat[k*dim3+i*dim2+j]  = U[i*dim2+j];
	  qmat_struct->v_mat[k*dim3+i*dim2+j]  = V[i*dim2+j];
	  if(V[i*dim2+j] > 1.E+20)
	    Exit("\n. Numerical overflow during the calculation of eigen vectors\n");
	}

      For(i,mod->ns)
	{
	  qmat_struct->expD_mr_vct[k*dim2+i] = (fit_double)exp(Root[i]);
	}
    }

  Free(U);
  Free(V);
  Free(Root);
}

/*********************************************************/

void Update_Eigen(int dim, fit_double *Q, fit_double *U, fit_double *V, fit_double *Root, model *mod)
{
    /* Call the eigenvalue decomposition routine with correction if
       round-off errors occured when inverting the instantaneous rate 
       matrix 
    */
    
  int i,j;
  int result;
  fit_double *di,*space,*Q_buff;
  fit_double scalar;
  int n_iter,*spaceint;
  fit_double **Q_buff_sq, **U_sq;
  fit_double *pi;
  
  di        = (fit_double *)mCalloc(dim,sizeof(fit_double));
  space     = (fit_double *)mCalloc(2*dim,sizeof(fit_double));
  spaceint  = (int    *)mCalloc(2*dim,sizeof(int));
  Q_buff    = (fit_double *)mCalloc(dim*dim,sizeof(fit_double));
  Q_buff_sq = (fit_double **)mCalloc(dim,sizeof(fit_double *));
  U_sq      = (fit_double **)mCalloc(dim,sizeof(fit_double *));
  pi        = (fit_double *)mCalloc(dim,sizeof(fit_double));
  
  For(i,dim) 
    {
      Q_buff_sq[i] = (fit_double *)mCalloc(dim,sizeof(fit_double ));
      U_sq[i]      = (fit_double *)mCalloc(dim,sizeof(fit_double ));
    }
  
  
  For(i,dim*dim) 
    {
      Q_buff[i] = Q[i];
      Q_buff_sq[i/dim][i%dim] = Q[i];
    }
  
  
  if(mod)
    {
      if((mod->model_number == 21) && (dim > 3))
	{
	  For(i,mod->n_omega)
	    For(j,mod->ns_codon)
	    {
	      pi[i*mod->ns_codon+j] =
		mod->pi[j] *
		mod->qmat_struct[0]->omega_proba[i];
	    }
	}
      else
	{
	  For(j,dim)
	    pi[j] = mod->pi[j];
	}
    }

  scalar   = 1.0;
  n_iter   = 0;
  result   = 0;
  
  if(!mod)
    {
      result = EigenRealGeneral(dim,Q_buff_sq,Root,di,U_sq,spaceint,space);
      For(i,dim*dim) U[i] = U_sq[i/dim][i%dim];
      /*       result = eigen(1,Q_buff,dim,Root,di,U,V,space); */
    }
  
  else
    {
      result = EigenRealGeneral(dim,Q_buff_sq,Root,di,U_sq,spaceint,space);
      For(i,dim*dim) U[i] = U_sq[i/dim][i%dim];
/*             result = eigenQREV(Q_buff,pi,mod->ns,Root,U,V,space); */
/*             result = eigen(1,Q_buff,dim,Root,di,U,V,space); */
    }

/*   For(i,dim*dim) printf("U[%d] = %E\n",i,U[i]);  */
/*   Exit(""); */

  if(!result)
    {      
      /* Compute inverse(Vr) into Vi */
      For (i,dim*dim) V[i] = U[i];
      while(!Matinv(V, dim, dim, space))
	{
	  printf("\n. Trying Q<-Q*scalar and then Root<-Root/scalar to fix this...\n");
	  fflush(NULL);
	  Q_buff   = (fit_double *)mCalloc(dim*dim,sizeof(fit_double));
	  scalar += scalar / 3.;
	  For(i,dim*dim) {Q_buff[i] = Q[i];}	  
	  For(i,dim*dim) Q_buff[i] *= scalar;
	  result = eigen(1,Q_buff,dim,Root,di,U,V,space);
	  if (result == -1)
	    Exit("\n. Eigenvalues/vectors computation did not converge : computation cancelled\n");
	  else if (result == 1)
	    Exit("\n. Complex eigenvalues/vectors : computation cancelled\n");
	  
	  For (i,dim*dim) V[i] = U[i];
	  n_iter++; 
	  if(n_iter > 100) Exit("\n. Cannot work out eigen vectors\n");
	};
      For(i,dim) Root[i] /= scalar;
    }
  
  Free(Q_buff);
  For(i,dim) { Free(Q_buff_sq[i]); Free(U_sq[i]);}
  Free(U_sq);
  Free(Q_buff_sq);
  Free(di);
  Free(space);
  Free(spaceint);
  Free(pi);
}

/*********************************************************/

void Init_T_Omega(model *mod)
{
  /*
    If switches (S1 or S2)

    n_catq
    ------->
 n |  o---|--->w0
 _ |______|
 o |      |
 m |  o------->w1
 e |______|
 g |      |
 a |  o------->w2
   |______|
   |
   V

   
   If no switches


                       n_omega
    ------------------------------------------------->                     
 n |  o---|--->w0	|  o---|--->NULL    |  o---|--->NULL    
 _ |______|		|______|            |______|          
 c |      |		|      |            |      |          
 a |  o------->w1	|  o------->NULL    |  o------->NULL  
 t |______|		|______|            |______|          
 q |      |		|      |            |      |          
   |  o------->w2	|  o------->NULL    |  o------->NULL  
   |______|          	|______|            |______|          
   |
   V

  */


  int i,j;

  For(i,mod->n_catq*mod->n_omega) mod->qmat_struct[0]->t_omega[i] = NULL;

  if(mod->switch_modelname == NO_SWITCH)
    {
      For(i,mod->n_catq)
	{
	  mod->qmat_struct[0]->t_omega[i*mod->n_omega + 0] = mod->qmat_struct[0]->omega+i;
	  for(j=1;j<mod->n_omega;j++)
	    mod->qmat_struct[0]->t_omega[i*mod->n_omega + j] = NULL;
	}
    }
  else
    {
      For(i,mod->n_catq)
	{
	  For(j,mod->n_omega)
	    mod->qmat_struct[0]->t_omega[i*mod->n_omega + j] = mod->qmat_struct[0]->omega+j;
	}
    }
}

/*********************************************************/

void Init_T_Omega_Proba(model *mod)
{
  /*
    If switches (S1 or S2)

    n_catq
    ------->
 n |  o---|--->p0
 _ |______|
 o |      |
 m |  o------->p1
 e |______|
 g |      |
 a |  o------->p2
   |______|
   |
   V

   
   If no switches


                       n_omega
    ------------------------------------------------->                     
 n |  o---|--->p0	|  o---|--->NULL    |  o---|--->NULL    
 _ |______|		|______|            |______|          
 c |      |		|      |            |      |          
 a |  o------->p1 	|  o------->NULL    |  o------->NULL  
 t |______|		|______|            |______|          
 q |      |		|      |            |      |          
   |  o------->p2	|  o------->NULL    |  o------->NULL  
   |______|          	|______|            |______|          
   |
   V

  */

  int i,j;


  if(mod->switch_modelname == NO_SWITCH)
    {
      For(i,mod->n_catq)
	{
	  mod->qmat_struct[0]->t_omega_proba[i*mod->n_omega + 0] = 
	    mod->qmat_struct[0]->omega_proba+i;
	  
	  for(j=1;j<mod->n_omega;j++)
	    mod->qmat_struct[0]->t_omega_proba[i*mod->n_omega + j] = NULL;
	}
    }
  else if((mod->switch_modelname == SWITCH_S1) || (mod->switch_modelname == SWITCH_S2))
    {
      For(i,mod->n_catq)
	For(j,mod->n_omega)
	mod->qmat_struct[0]->t_omega_proba[i*mod->n_omega + j] = 
	mod->qmat_struct[0]->omega_proba+j;
    }
  else
    Exit("Try again...\n");
}

/*********************************************************/
